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