Public Member Functions | Private Member Functions | Private Attributes | List of all members
ctp::CTPHelper Class Reference

Class containing some utility functions for all things CVN. More...

#include <CTPHelper.h>

Public Member Functions

 CTPHelper (const fhicl::ParameterSet &pset)
 
 ~CTPHelper ()
 
const ctp::CTPResult RunConvolutionalTrackPID (const art::Ptr< recob::PFParticle > particle, const art::Event &evt) const
 
const std::vector< std::vector< float > > GetNetworkInputs (const art::Ptr< recob::PFParticle >, const art::Event &evt) const
 
const std::vector< float > GetDeDxVector (const art::Ptr< recob::PFParticle >, const art::Event &evt) const
 
const std::vector< float > GetVariableVector (const art::Ptr< recob::PFParticle >, const art::Event &evt) const
 
const std::pair< const simb::MCParticle *, float > GetTrueParticle (const art::Ptr< recob::PFParticle >, const art::Event &evt) const
 
const int GetTruePDGCode (const art::Ptr< recob::PFParticle >, const art::Event &evt) const
 

Private Member Functions

void SmoothDedxVector (std::vector< float > &dedx) const
 
void PadDedxVector (std::vector< float > &dedx, const float mean, const float sigma) const
 
void GetDedxMeanAndSigma (const std::vector< float > &dedx, float &mean, float &sigma) const
 
void GetDeflectionMeanAndSigma (const art::Ptr< recob::Track > track, float &mean, float &sigma) const
 
void GetChildParticles (const art::Ptr< recob::PFParticle > part, const art::Event &evt, float &nTrack, float &nShower, float &nGrand) const
 
void NormaliseInputs (std::vector< std::vector< float >> &netInputs) const
 

Private Attributes

std::string fNetDir
 
std::string fNetName
 
std::string fParticleLabel
 
std::string fTrackLabel
 
std::string fShowerLabel
 
std::string fCalorimetryLabel
 
unsigned int fMinTrackPoints
 
unsigned int fDedxLength
 
float fQMax
 
float fQJump
 
bool fNormalise
 

Detailed Description

Class containing some utility functions for all things CVN.

Definition at line 30 of file CTPHelper.h.

Constructor & Destructor Documentation

ctp::CTPHelper::CTPHelper ( const fhicl::ParameterSet pset)

Definition at line 39 of file CTPHelper.cxx.

39  {
40  fNetDir = pset.get<std::string>("NetworkPath");
41  fNetName = pset.get<std::string>("NetworkName");
42  fParticleLabel = pset.get<std::string>("ParticleLabel");
43  fTrackLabel = pset.get<std::string>("TrackLabel");
44  fShowerLabel = pset.get<std::string>("ShowerLabel");
45  fCalorimetryLabel = pset.get<std::string>("CalorimetryLabel");
46  fMinTrackPoints = pset.get<unsigned int>("MinHits",50);
47  fDedxLength = pset.get<unsigned int>("DedxLength",100);
48  fQMax = pset.get<float>("MaxCharge",1000);
49  fQJump = pset.get<float>("MaxChargeJump",500);
50  fNormalise = pset.get<bool>("NormaliseInputs",true);
51  }
std::string fShowerLabel
Definition: CTPHelper.h:66
unsigned int fDedxLength
Definition: CTPHelper.h:71
std::string fNetDir
Definition: CTPHelper.h:60
bool fNormalise
Definition: CTPHelper.h:75
std::string string
Definition: nybbler.cc:12
std::string fTrackLabel
Definition: CTPHelper.h:65
T get(std::string const &key) const
Definition: ParameterSet.h:271
unsigned int fMinTrackPoints
Definition: CTPHelper.h:70
std::string fNetName
Definition: CTPHelper.h:61
float fQJump
Definition: CTPHelper.h:73
std::string fParticleLabel
Definition: CTPHelper.h:64
std::string fCalorimetryLabel
Definition: CTPHelper.h:67
ctp::CTPHelper::~CTPHelper ( )

Definition at line 53 of file CTPHelper.cxx.

53  {
54  }

Member Function Documentation

void ctp::CTPHelper::GetChildParticles ( const art::Ptr< recob::PFParticle part,
const art::Event evt,
float &  nTrack,
float &  nShower,
float &  nGrand 
) const
private

Definition at line 276 of file CTPHelper.cxx.

276  {
277 
278  nTrack = 0.;
279  nShower = 0.;
280  nGrand = 0.;
281 
282  std::vector<art::Ptr<recob::PFParticle>> children = dune_ana::DUNEAnaPFParticleUtils::GetChildParticles(part,evt,fParticleLabel);
283 
284  for(const art::Ptr<recob::PFParticle> child : children){
287  nGrand += child->NumDaughters();
288  }
289 // std::cout << "Children = " << children.size() << "( " << nTrack << ", " << nShower << ") and grand children = " << nGrand << std::endl;
290  }
std::string fShowerLabel
Definition: CTPHelper.h:66
static bool IsShower(const art::Ptr< recob::PFParticle > &pParticle, const art::Event &evt, const std::string &pParticleLabel, const std::string &showerLabel)
Check if this particle has an associated shower.
std::string fTrackLabel
Definition: CTPHelper.h:65
static std::vector< art::Ptr< recob::PFParticle > > GetChildParticles(const art::Ptr< recob::PFParticle > &pParticle, const art::Event &evt, const std::string &label)
Get the child particles (one step down in the hierarchy) of this particle.
std::string fParticleLabel
Definition: CTPHelper.h:64
static bool IsTrack(const art::Ptr< recob::PFParticle > &pParticle, const art::Event &evt, const std::string &pParticleLabel, const std::string &trackLabel)
Check if this particle has an associated track.
void ctp::CTPHelper::GetDedxMeanAndSigma ( const std::vector< float > &  dedx,
float &  mean,
float &  sigma 
) const
private

Definition at line 246 of file CTPHelper.cxx.

246  {
247 
248  float averageDedx = 0;
249  float sigmaDedx = 0;
250  for(const float &q : dedx) averageDedx += q;
251  mean = averageDedx / static_cast<float>(dedx.size());
252  for(const float &q : dedx) sigmaDedx += (mean -q)*(mean-q);
253  sigma = std::sqrt(sigmaDedx / static_cast<float>(dedx.size()));
254  }
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:16
const std::vector< float > ctp::CTPHelper::GetDeDxVector ( const art::Ptr< recob::PFParticle part,
const art::Event evt 
) const

Definition at line 144 of file CTPHelper.cxx.

144  {
145  return GetNetworkInputs(part,evt).at(0);
146  }
const std::vector< std::vector< float > > GetNetworkInputs(const art::Ptr< recob::PFParticle >, const art::Event &evt) const
Definition: CTPHelper.cxx:81
void ctp::CTPHelper::GetDeflectionMeanAndSigma ( const art::Ptr< recob::Track track,
float &  mean,
float &  sigma 
) const
private

Definition at line 256 of file CTPHelper.cxx.

256  {
257 
258  std::vector<float> trajAngle;
259  for(unsigned int p = 1; p < track->Trajectory().NPoints(); ++p){
260  TVector3 thisDir = track->Trajectory().DirectionAtPoint<TVector3>(p);
261  TVector3 prevDir = track->Trajectory().DirectionAtPoint<TVector3>(p-1);
262  trajAngle.push_back(thisDir.Angle(prevDir));
263 
264 // if(p < 11 || p > track->Trajectory().NPoints() - 11) std::cout << "Deflection at point " << p << " = " << trajAngle.at(p-1) << std::endl;
265  }
266 
267  // Average and sigma of the angular deflection between trajectory points (wobble)
268  float averageAngle = 0;
269  float standardDevAngle = 0;
270  for(const float &a : trajAngle) averageAngle += a;
271  mean = averageAngle / static_cast<float>(trajAngle.size());
272  for(const float &a : trajAngle) standardDevAngle += (a - mean)*(a - mean);
273  sigma = sqrt(standardDevAngle / static_cast<float>(trajAngle.size()));
274  }
T DirectionAtPoint(unsigned int p) const
Direction at point p. Use e.g. as:
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
Definition: Track.h:98
size_t NPoints() const
Returns the number of stored trajectory points.
Definition: Trajectory.h:167
const double a
p
Definition: test.py:223
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:16
const std::vector< std::vector< float > > ctp::CTPHelper::GetNetworkInputs ( const art::Ptr< recob::PFParticle part,
const art::Event evt 
) const

Definition at line 81 of file CTPHelper.cxx.

81  {
82 
83  std::vector<std::vector<float>> netInputs;
84 
86 // std::cout << "CTPHelper: this PFParticle is not track-like... returning empty vector." << std::endl;
87  return std::vector<std::vector<float>>();
88  }
89 
90  // Use the analysis utilities to simplify finding products and associations
93 
94  if(thisCalo->dEdx().size() < fMinTrackPoints){
95 // std::cout << "CTPHelper: this track has too few points for PID (" << thisCalo->dEdx().size() << ")... returning empty vector." << std::endl;
96  return std::vector<std::vector<float>>();
97  }
98 
99  std::vector<float> dedxVector = thisCalo->dEdx();
100  this->SmoothDedxVector(dedxVector);
101  float dedxMean = 0.;
102  float dedxSigma = 0.;
103 
104  // We want to use the middle third of the dedx vector (with max length 100)
105  std::vector<float> dedxTrunc;
106  unsigned int pointsForAverage = (fDedxLength - fMinTrackPoints) / 3;
107  unsigned int avStart = dedxVector.size() - 1 - pointsForAverage;
108  unsigned int avEnd = dedxVector.size() - 1 - (2*pointsForAverage);
109  for(unsigned int e = avStart; e > avEnd; --e) dedxTrunc.push_back(dedxVector.at(e));
110  this->GetDedxMeanAndSigma(dedxTrunc,dedxMean,dedxSigma);
111 
112  // If our dedx vector is between fMinTrackPoints and fDedxLength in size then we need to pad it
113  if(dedxVector.size() < fDedxLength){
114  this->PadDedxVector(dedxVector,dedxMean,dedxSigma);
115  }
116 
117  std::vector<float> finalInputDedx;
118  finalInputDedx.insert(finalInputDedx.begin(),dedxVector.end() - fDedxLength,dedxVector.end());
119 
120  std::vector<float> finalInputVariables;
121  // Get the number of child particles
122  float nTrack, nShower, nGrand;
123  this->GetChildParticles(part,evt,nTrack,nShower,nGrand);
124  finalInputVariables.push_back(nTrack);
125  finalInputVariables.push_back(nShower);
126  finalInputVariables.push_back(nGrand);
127  // Now add the dedx mean and sigma
128  finalInputVariables.push_back(dedxMean);
129  finalInputVariables.push_back(dedxSigma);
130  // Finally, get the angular deflection mean and sigma
131  float deflectionMean, deflectionSigma;
132  this->GetDeflectionMeanAndSigma(thisTrack,deflectionMean,deflectionSigma);
133  finalInputVariables.push_back(deflectionMean);
134  finalInputVariables.push_back(deflectionSigma);
135 
136  netInputs.push_back(finalInputDedx);
137  netInputs.push_back(finalInputVariables);
138 
139  if(fNormalise) this->NormaliseInputs(netInputs);
140 
141  return netInputs;
142  }
unsigned int fDedxLength
Definition: CTPHelper.h:71
static art::Ptr< recob::Track > GetTrack(const art::Ptr< recob::PFParticle > &pParticle, const art::Event &evt, const std::string &pParticleLabel, const std::string &trackLabel)
Get the track associated to this particle. Should only be called if IsTrack method succeeds...
void NormaliseInputs(std::vector< std::vector< float >> &netInputs) const
Definition: CTPHelper.cxx:292
bool fNormalise
Definition: CTPHelper.h:75
std::string fTrackLabel
Definition: CTPHelper.h:65
void SmoothDedxVector(std::vector< float > &dedx) const
Definition: CTPHelper.cxx:189
const double e
void GetDeflectionMeanAndSigma(const art::Ptr< recob::Track > track, float &mean, float &sigma) const
Definition: CTPHelper.cxx:256
void GetDedxMeanAndSigma(const std::vector< float > &dedx, float &mean, float &sigma) const
Definition: CTPHelper.cxx:246
const std::vector< float > & dEdx() const
Definition: Calorimetry.h:101
unsigned int fMinTrackPoints
Definition: CTPHelper.h:70
std::string fParticleLabel
Definition: CTPHelper.h:64
static art::Ptr< anab::Calorimetry > GetCalorimetry(const art::Ptr< recob::Track > &pTrack, const art::Event &evt, const std::string &trackLabel, const std::string &caloLabel)
Get the particle associated with the track.
static bool IsTrack(const art::Ptr< recob::PFParticle > &pParticle, const art::Event &evt, const std::string &pParticleLabel, const std::string &trackLabel)
Check if this particle has an associated track.
std::string fCalorimetryLabel
Definition: CTPHelper.h:67
void GetChildParticles(const art::Ptr< recob::PFParticle > part, const art::Event &evt, float &nTrack, float &nShower, float &nGrand) const
Definition: CTPHelper.cxx:276
void PadDedxVector(std::vector< float > &dedx, const float mean, const float sigma) const
Definition: CTPHelper.cxx:224
const std::pair< const simb::MCParticle *, float > ctp::CTPHelper::GetTrueParticle ( const art::Ptr< recob::PFParticle part,
const art::Event evt 
) const

Definition at line 152 of file CTPHelper.cxx.

152  {
153  // Get hits
154  const std::vector<art::Ptr<recob::Hit>> collectionHits = dune_ana::DUNEAnaPFParticleUtils::GetHits(part,evt,fParticleLabel);
155  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
156 
157  // Function to find a weighted and sorted vector of the true particles matched to a particle
158  using weightedMCPair = std::pair<const simb::MCParticle*, float>;
159  std::vector<weightedMCPair> outVecHits;
160 
161  // Loop over all hits in the input vector and record the contributing MCParticles.
164  std::unordered_map<const simb::MCParticle*, float> mcHitMap;
165  float hitTotal = 0;
166  for(const art::Ptr<recob::Hit> &hit : collectionHits) {
167  for(const sim::TrackIDE& ide : bt_serv->HitToTrackIDEs(clockData, hit)) {
168  const simb::MCParticle* curr_part = pi_serv->TrackIdToParticle_P(ide.trackID);
169  mcHitMap[curr_part] += 1.;
170  ++hitTotal;
171  }
172  }
173 
174  for (weightedMCPair const &p : mcHitMap) outVecHits.push_back(p);
175  // Can't continue without a truth match
176  if(outVecHits.size() == 0) assert(0);
177 
178  std::sort(outVecHits.begin(),outVecHits.end(),[](weightedMCPair a, weightedMCPair b){ return a.second > b.second;});
179  for(weightedMCPair &p : outVecHits) p.second /= hitTotal;
180 
181 // std::cout << "Returning truth match..." << std::endl;
182  return outVecHits.at(0);
183  }
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
const simb::MCParticle * TrackIdToParticle_P(int id) const
const double a
p
Definition: test.py:223
Detector simulation of raw signals on wires.
std::string fParticleLabel
Definition: CTPHelper.h:64
static std::vector< art::Ptr< recob::Hit > > GetHits(const art::Ptr< recob::PFParticle > &pParticle, const art::Event &evt, const std::string &label)
Get the hits associated to this particle.
static bool * b
Definition: config.cpp:1043
Ionization energy from a Geant4 track.
Definition: SimChannel.h:25
const int ctp::CTPHelper::GetTruePDGCode ( const art::Ptr< recob::PFParticle part,
const art::Event evt 
) const

Definition at line 185 of file CTPHelper.cxx.

185  {
186  return this->GetTrueParticle(part,evt).first->PdgCode();
187  }
const std::pair< const simb::MCParticle *, float > GetTrueParticle(const art::Ptr< recob::PFParticle >, const art::Event &evt) const
Definition: CTPHelper.cxx:152
const std::vector< float > ctp::CTPHelper::GetVariableVector ( const art::Ptr< recob::PFParticle part,
const art::Event evt 
) const

Definition at line 148 of file CTPHelper.cxx.

148  {
149  return GetNetworkInputs(part,evt).at(1);
150  }
const std::vector< std::vector< float > > GetNetworkInputs(const art::Ptr< recob::PFParticle >, const art::Event &evt) const
Definition: CTPHelper.cxx:81
void ctp::CTPHelper::NormaliseInputs ( std::vector< std::vector< float >> &  netInputs) const
private

Definition at line 292 of file CTPHelper.cxx.

292  {
293 
294  // Firstly, we need to normalise the dE/dx values
295  for(float &dedx : inputs.at(0)){
296  // Protect against negative values
297  if(dedx < 1.e-3){
298 // std::cout << "Very low dedx value = " << dedx << std::endl;
299  dedx = 1.e-3;
300  }
301  float newVal = std::log(2*dedx) + 1;
302  if(newVal > 5.) newVal = 5.;
303  dedx = (newVal / 2.5) - 1;
304  }
305 
306  // For the three child values
307  for(unsigned int v = 0; v < 3; ++v){
308  float val = inputs.at(1).at(v);
309  if(val > 5) val = 5;
310  inputs.at(1).at(v) = (val / 2.5) - 1.0;
311  }
312 
313  // The charge mean
314  float val = inputs.at(1).at(3);
315  val = std::log(2*val) + 1;
316  if(val > 5.) val = 5.;
317  inputs.at(1).at(3) = (val / 2.5) - 1;
318 
319  // The charge sigma
320  val = inputs.at(1).at(4);
321  if(val > 3) val = 3.;
322  inputs.at(1).at(4) = (val/1.5) - 1;
323 
324  // The angle mean
325  val = inputs.at(1).at(5);
326  if(val > 0.05) val = 0.05;
327  inputs.at(1).at(5) = (val / 0.025) - 1.0;
328 
329  // The angle sigma
330  val = inputs.at(1).at(6);
331  if(val > 0.3) val = 0.3;
332  inputs.at(1).at(6) = (val / 0.15) - 1.0;
333 
334  }
const double e
void ctp::CTPHelper::PadDedxVector ( std::vector< float > &  dedx,
const float  mean,
const float  sigma 
) const
private

Definition at line 224 of file CTPHelper.cxx.

224  {
225 
226  std::default_random_engine generator;
227  std::normal_distribution<float> gaussDist(mean,sigma);
228 
229  unsigned int originalSize = dedx.size();
230 
231  for (unsigned int h = 0; h + originalSize < fDedxLength; ++h)
232  {
233  // Pick a random Gaussian value but ensure we don't go negative
234  float randVal = -1;
235  do
236  {
237  randVal = gaussDist(generator);
238  }
239  while (randVal < 0);
240  // Pad from beginning to keep the real track part at the end
241  dedx.insert(dedx.begin(),randVal);
242  }
243 
244  }
unsigned int fDedxLength
Definition: CTPHelper.h:71
generator
Definition: train.py:468
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:16
const ctp::CTPResult ctp::CTPHelper::RunConvolutionalTrackPID ( const art::Ptr< recob::PFParticle particle,
const art::Event evt 
) const

Definition at line 57 of file CTPHelper.cxx.

57  {
58 
59  // Get the inputs to the network
60  std::vector< std::vector< std::vector<float> > > finalInputs;
61 
62  std::vector< std::vector<float> > twoVecs = GetNetworkInputs(part,evt);
63  if(twoVecs.empty()) return ctp::CTPResult();
64 
65  finalInputs.push_back(twoVecs);
66 
67 // std::cout << "We have sensible inputs... building TF object" << std::endl;
68  // Load the network and run it
69  const std::string fullPath = cet::getenv(fNetDir) + "/" + fNetName;
70  std::unique_ptr<tf::CTPGraph> convNet = tf::CTPGraph::create(fullPath.c_str(),std::vector<std::string>(),2,1);
71 // std::cout << "Calling TF interface" << std::endl;
72  std::vector< std::vector< std::vector<float> > > convNetOutput = convNet->run(finalInputs);
73 
74 // std::cout << "Got result from TF" << std::endl;
75  ctp::CTPResult result(convNetOutput.at(0).at(0));
76 // std::cout << "Returning CTPResult" << std::endl;
77  return result;
78  }
std::string fNetDir
Definition: CTPHelper.h:60
static QCString result
Class containing some utility functions for all things CVN.
Definition: CTPResult.h:16
static std::unique_ptr< CTPGraph > create(const char *graph_file_name, const std::vector< std::string > &outputs={}, int ninputs=1, int noutputs=1)
Definition: CTPGraph.h:32
std::string string
Definition: nybbler.cc:12
const std::vector< std::vector< float > > GetNetworkInputs(const art::Ptr< recob::PFParticle >, const art::Event &evt) const
Definition: CTPHelper.cxx:81
std::vector< float > run(const std::vector< std::vector< float > > &x)
std::string getenv(std::string const &name)
Definition: getenv.cc:15
std::string fNetName
Definition: CTPHelper.h:61
void ctp::CTPHelper::SmoothDedxVector ( std::vector< float > &  dedx) const
private

Definition at line 189 of file CTPHelper.cxx.

189  {
190 
191  // Firstly, get rid of all very high values > fQMax
192  for(float &val : dedx){
193 // std::cout << "Checking de/dx = " << val;
194  if(val > fQMax){
195 // std::cout << "CTPHelper: large value " << val << " set to " << fQMax << std::endl;
196  val = fQMax;
197  }
198  if(val < 1e-3){
199  std::cout << "CTPHelper: small value " << val << " set to 1.0e-3 " << std::endl;
200  val = 1e-3;
201  }
202 // std::cout << " value becomes " << val << std::endl;
203  }
204 
205  // Now try to smooth over jumps
206  unsigned int nQ = dedx.size();
207  // Now do the rest of the points
208  for(unsigned int q = 1; q < nQ - 1; ++q){
209  if((dedx[q] - dedx[q-1]) > fQJump)
210  {
211 // std::cout << "Updating dedx point " << q << " from " << dedx[q] << " to " << 0.5 * (dedx[q-1]+dedx[q+1]) << std::endl;
212  dedx[q] = 0.5 * (dedx[q-1]+dedx[q+1]);
213  }
214  }
215  // First and last points are special cases
216 // std::cout << "Initial first and last points: " << dedx[0] << ", " << dedx[nQ-1] << std::endl;
217 // if((dedx[0] - dedx[1]) > fQJump) dedx[0] = dedx[1] + (dedx[1] - dedx[2]);
218 // if((dedx[nQ-1] - dedx[nQ-2]) > fQJump) dedx[nQ-1] = dedx[nQ-2] + (dedx[nQ-2] - dedx[nQ-3]);
219 // std::cout << "Final first and last points: " << dedx[0] << ", " << dedx[nQ-1] << std::endl;
220 // std::cout << "Using... " << dedx[1] << ", " << dedx[nQ-3] << ", " << dedx[nQ-2] << std::endl;
221 
222  }
const double e
float fQJump
Definition: CTPHelper.h:73
QTextStream & endl(QTextStream &s)

Member Data Documentation

std::string ctp::CTPHelper::fCalorimetryLabel
private

Definition at line 67 of file CTPHelper.h.

unsigned int ctp::CTPHelper::fDedxLength
private

Definition at line 71 of file CTPHelper.h.

unsigned int ctp::CTPHelper::fMinTrackPoints
private

Definition at line 70 of file CTPHelper.h.

std::string ctp::CTPHelper::fNetDir
private

Definition at line 60 of file CTPHelper.h.

std::string ctp::CTPHelper::fNetName
private

Definition at line 61 of file CTPHelper.h.

bool ctp::CTPHelper::fNormalise
private

Definition at line 75 of file CTPHelper.h.

std::string ctp::CTPHelper::fParticleLabel
private

Definition at line 64 of file CTPHelper.h.

float ctp::CTPHelper::fQJump
private

Definition at line 73 of file CTPHelper.h.

float ctp::CTPHelper::fQMax
private

Definition at line 72 of file CTPHelper.h.

std::string ctp::CTPHelper::fShowerLabel
private

Definition at line 66 of file CTPHelper.h.

std::string ctp::CTPHelper::fTrackLabel
private

Definition at line 65 of file CTPHelper.h.


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