Classes | Functions
pduneana Namespace Reference

Classes

struct  calo_point
 
struct  cnnOutput2D
 
class  PDSPAnalyzer
 

Functions

double distance2 (double x, double y, double z, double *p)
 
void line (double t, double *p, double &x, double &y, double &z)
 
void SumDistance2 (int &, double *, double &sum, double *par, int)
 
TVector3 FitLine (const std::vector< TVector3 > &input)
 
bool sort_IDEs (const sim::IDE *i1, const sim::IDE *i2)
 
std::map< int, std::vector< const sim::IDE * > > slice_IDEs (std::vector< const sim::IDE * > ides, double the_z0, double the_pitch, double true_endZ)
 
TFile * OpenFile (const std::string filename)
 
cnnOutput2D GetCNNOutputFromPFParticle (const recob::PFParticle &part, const art::Event &evt, const anab::MVAReader< recob::Hit, 4 > &CNN_results, protoana::ProtoDUNEPFParticleUtils &pfpUtil, std::string fPFParticleTag)
 
cnnOutput2D GetCNNOutputFromPFParticleFromPlane (const recob::PFParticle &part, const art::Event &evt, const anab::MVAReader< recob::Hit, 4 > &CNN_results, protoana::ProtoDUNEPFParticleUtils &pfpUtil, std::string fPFParticleTag, size_t planeID)
 

Function Documentation

double pduneana::distance2 ( double  x,
double  y,
double  z,
double *  p 
)

Definition at line 4752 of file PDSPAnalyzer_module.cc.

4752  {
4753  // distance line point is D= | (xp-x0) cross ux |
4754  // where ux is direction of line and x0 is a point in the line (like t = 0)
4755  ROOT::Math::XYZVector xp(x,y,z);
4756  ROOT::Math::XYZVector x0(p[0], p[2], 0.);
4757  ROOT::Math::XYZVector x1(p[0] + p[1], p[2] + p[3], 1.);
4758  ROOT::Math::XYZVector u = (x1-x0).Unit();
4759  double d2 = ((xp-x0).Cross(u)) .Mag2();
4760  return d2;
4761 }
p
Definition: test.py:223
list x
Definition: train.py:276
TVector3 pduneana::FitLine ( const std::vector< TVector3 > &  input)

Definition at line 4780 of file PDSPAnalyzer_module.cc.

4780  {
4781  TGraph2D * gr = new TGraph2D();
4782  for (size_t i = 0; i < input.size(); ++i) {
4783  gr->SetPoint(i, input[i].X(), input[i].Y(), input[i].Z());
4784  }
4785 
4786  TVirtualFitter * min = TVirtualFitter::Fitter(0,4);
4787  min->SetObjectFit(gr);
4788  min->SetFCN(SumDistance2);
4789 
4790  double arglist[10];
4791 
4792  arglist[0] = -1;
4793  min->ExecuteCommand("SET PRINT",arglist,1);
4794 
4795 
4796  double pStart[4] = {1,1,1,1};
4797  min->SetParameter(0,"x0",pStart[0],0.01,0,0);
4798  min->SetParameter(1,"Ax",pStart[1],0.01,0,0);
4799  min->SetParameter(2,"y0",pStart[2],0.01,0,0);
4800  min->SetParameter(3,"Ay",pStart[3],0.01,0,0);
4801 
4802  arglist[0] = 1000; // number of function calls
4803  arglist[1] = 0.001; // tolerance
4804  min->ExecuteCommand("MIGRAD", arglist, 2);
4805 
4806  // get fit parameters
4807  double parFit[4];
4808  for (int i = 0; i < 4; ++i) {
4809  parFit[i] = min->GetParameter(i);
4810  }
4811  double startX1, startY1, startZ1;
4812  double startX2, startY2, startZ2;
4813  line(0, parFit, startX1, startY1, startZ1);
4814  line(1, parFit, startX2, startY2, startZ2);
4815 
4816  TVector3 diff(startX2 - startX1,
4817  startY2 - startY1,
4818  startZ2 - startZ1);
4819  delete gr;
4820  delete min;
4821  return diff;
4822 }
void SumDistance2(int &, double *, double &sum, double *par, int)
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
void line(double t, double *p, double &x, double &y, double &z)
cnnOutput2D pduneana::GetCNNOutputFromPFParticle ( const recob::PFParticle part,
const art::Event evt,
const anab::MVAReader< recob::Hit, 4 > &  CNN_results,
protoana::ProtoDUNEPFParticleUtils pfpUtil,
std::string  fPFParticleTag 
)

Definition at line 197 of file PDSPAnalyzer_module.cc.

201  {
202 
203  cnnOutput2D output;
204  const std::vector< art::Ptr< recob::Hit > > daughterPFP_hits
205  = pfpUtil.GetPFParticleHits_Ptrs(part, evt, fPFParticleTag);
206 
207  double tot_charge = 0.;
208  for (size_t h = 0; h < daughterPFP_hits.size(); ++h) {
209  std::array<float,4> cnn_out = CNN_results.getOutput( daughterPFP_hits[h] );
210  double hitcharge = daughterPFP_hits[h]->Integral();
211  double track_score = cnn_out[ CNN_results.getIndex("track") ];
212  double em_score = cnn_out[ CNN_results.getIndex("em") ];
213  double michel_score = cnn_out[ CNN_results.getIndex("michel") ];
214  double none_score = cnn_out[ CNN_results.getIndex("none") ];
215  output.track += track_score;
216  output.em += em_score;
217  output.michel += michel_score;
218  output.none += none_score;
219  output.track_weight_by_charge += hitcharge*track_score;
220  output.em_weight_by_charge += hitcharge*em_score;
221  output.michel_weight_by_charge += hitcharge*michel_score;
222  output.none_weight_by_charge += hitcharge*none_score;
223  tot_charge += hitcharge;
224  }
225  output.nHits = daughterPFP_hits.size();
226  if (tot_charge != 0) {
227  output.track_weight_by_charge /= tot_charge;
228  output.em_weight_by_charge /= tot_charge;
229  output.michel_weight_by_charge /= tot_charge;
230  output.none_weight_by_charge /= tot_charge;
231  }
232  else {
233  output.track_weight_by_charge = -999.;
234  output.em_weight_by_charge = -999.;
235  output.michel_weight_by_charge = -999.;
236  output.none_weight_by_charge = -999.;
237  }
238 
239  return output;
240  }
int getIndex(const std::string &name) const
Index of column with given name, or -1 if name not found.
Definition: MVAReader.h:82
const std::vector< art::Ptr< recob::Hit > > GetPFParticleHits_Ptrs(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel) const
std::array< float, N > getOutput(size_t key) const
Get copy of the MVA output vector at index "key".
Definition: MVAReader.h:129
cnnOutput2D pduneana::GetCNNOutputFromPFParticleFromPlane ( const recob::PFParticle part,
const art::Event evt,
const anab::MVAReader< recob::Hit, 4 > &  CNN_results,
protoana::ProtoDUNEPFParticleUtils pfpUtil,
std::string  fPFParticleTag,
size_t  planeID 
)

Definition at line 244 of file PDSPAnalyzer_module.cc.

244  {
245 
246  cnnOutput2D output;
247  const std::vector< art::Ptr< recob::Hit > > daughterPFP_hits = pfpUtil.GetPFParticleHitsFromPlane_Ptrs( part, evt, fPFParticleTag, planeID );
248 
249  double tot_charge = 0.;
250  for (size_t h = 0; h < daughterPFP_hits.size(); ++h) {
251  std::array<float,4> cnn_out = CNN_results.getOutput( daughterPFP_hits[h] );
252  double hitcharge = daughterPFP_hits[h]->Integral();
253  double track_score = cnn_out[ CNN_results.getIndex("track") ];
254  double em_score = cnn_out[ CNN_results.getIndex("em") ];
255  double michel_score = cnn_out[ CNN_results.getIndex("michel") ];
256  double none_score = cnn_out[ CNN_results.getIndex("none") ];
257  output.track += track_score;
258  output.em += em_score;
259  output.michel += michel_score;
260  output.none += none_score;
261  output.track_weight_by_charge += hitcharge*track_score;
262  output.em_weight_by_charge += hitcharge*em_score;
263  output.michel_weight_by_charge += hitcharge*michel_score;
264  output.none_weight_by_charge += hitcharge*none_score;
265  tot_charge += hitcharge;
266  }
267  output.nHits = daughterPFP_hits.size();
268  if (tot_charge != 0) {
269  output.track_weight_by_charge /= tot_charge;
270  output.em_weight_by_charge /= tot_charge;
271  output.michel_weight_by_charge /= tot_charge;
272  output.none_weight_by_charge /= tot_charge;
273  }
274  else {
275  output.track_weight_by_charge = -999.;
276  output.em_weight_by_charge = -999.;
277  output.michel_weight_by_charge = -999.;
278  output.none_weight_by_charge = -999.;
279  }
280 
281  return output;
282  }
int getIndex(const std::string &name) const
Index of column with given name, or -1 if name not found.
Definition: MVAReader.h:82
const std::vector< art::Ptr< recob::Hit > > GetPFParticleHitsFromPlane_Ptrs(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, size_t planeID) const
std::array< float, N > getOutput(size_t key) const
Get copy of the MVA output vector at index "key".
Definition: MVAReader.h:129
void pduneana::line ( double  t,
double *  p,
double &  x,
double &  y,
double &  z 
)

Definition at line 4741 of file PDSPAnalyzer_module.cc.

4742  {
4743  // a parameteric line is define from 6 parameters but 4 are independent
4744  // x0,y0,z0,z1,y1,z1 which are the coordinates of two points on the line
4745  // can choose z0 = 0 if line not parallel to x-y plane and z1 = 1;
4746  x = p[0] + p[1]*t;
4747  y = p[2] + p[3]*t;
4748  z = t;
4749 }
p
Definition: test.py:223
list x
Definition: train.py:276
TFile* pduneana::OpenFile ( const std::string  filename)

Definition at line 121 of file PDSPAnalyzer_module.cc.

121  {
122  TFile * theFile = 0x0;
123  mf::LogInfo("pduneana::OpenFile") << "Searching for " << filename;
124  if (cet::file_exists(filename)) {
125  mf::LogInfo("pduneana::OpenFile") << "File exists. Opening " << filename;
126  theFile = new TFile(filename.c_str());
127  if (!theFile ||theFile->IsZombie() || !theFile->IsOpen()) {
128  delete theFile;
129  theFile = 0x0;
130  throw cet::exception("PDSPAnalyzer_module.cc") << "Could not open " << filename;
131  }
132  }
133  else {
134  mf::LogInfo("pduneana::OpenFile") << "File does not exist here. Searching FW_SEARCH_PATH";
135  cet::search_path sp{"FW_SEARCH_PATH"};
136  std::string found_filename;
137  auto found = sp.find_file(filename, found_filename);
138  if (!found) {
139  throw cet::exception("PDSPAnalyzer_module.cc") << "Could not find " << filename;
140  }
141 
142  mf::LogInfo("pduneana::OpenFile") << "Found file " << found_filename;
143  theFile = new TFile(found_filename.c_str());
144  if (!theFile ||theFile->IsZombie() || !theFile->IsOpen()) {
145  delete theFile;
146  theFile = 0x0;
147  throw cet::exception("PDSPAnalyzer_module.cc") << "Could not open " << found_filename;
148  }
149  }
150  return theFile;
151  }
std::string string
Definition: nybbler.cc:12
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
string filename
Definition: train.py:213
bool file_exists(std::string const &qualified_filename)
Definition: filesystem.cc:14
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::map<int, std::vector<const sim::IDE*> > pduneana::slice_IDEs ( std::vector< const sim::IDE * >  ides,
double  the_z0,
double  the_pitch,
double  true_endZ 
)

Definition at line 99 of file PDSPAnalyzer_module.cc.

101  {
102 
103  std::map<int, std::vector<const sim::IDE*>> results;
104 
105  for (size_t i = 0; i < ides.size(); ++i) {
106  int slice_num = std::floor(
107  (ides[i]->z - (the_z0 - the_pitch/2.)) / the_pitch);
108  /*
109  std::cout << "IDE: " << i << " ID: " << ides[i]->trackID << " Edep: "
110  << ides[i]->energy << " (X,Y,Z): " << "(" << ides[i]->x << ","
111  << ides[i]->y<<","<<ides[i]->z << ") Z0: " << the_z0
112  << " Slice: " << slice_num << std::endl;*/
113  results[slice_num].push_back(ides[i]);
114  }
115 
116  return results;
117  }
bool pduneana::sort_IDEs ( const sim::IDE i1,
const sim::IDE i2 
)

Definition at line 96 of file PDSPAnalyzer_module.cc.

96  {
97  return( i1->z < i2->z );
98  }
float z
z position of ionization [cm]
Definition: SimChannel.h:117
void pduneana::SumDistance2 ( int &  ,
double *  ,
double &  sum,
double *  par,
int   
)

Definition at line 4764 of file PDSPAnalyzer_module.cc.

4765  {
4766  // the TGraph must be a global variable
4767  TGraph2D * gr = dynamic_cast<TGraph2D*>( (TVirtualFitter::GetFitter())->GetObjectFit() );
4768  assert(gr != 0);
4769  double * x = gr->GetX();
4770  double * y = gr->GetY();
4771  double * z = gr->GetZ();
4772  int npoints = gr->GetN();
4773  sum = 0;
4774  for (int i = 0; i < npoints; ++i) {
4775  double d = distance2(x[i],y[i],z[i],par);
4776  sum += d;
4777  }
4778 }
list x
Definition: train.py:276
double distance2(double x, double y, double z, double *p)