TrackCleanUpAlg.cxx
Go to the documentation of this file.
1 
2 ////////////////////////////////////////////////////////////////////////
3 /// \file TrackCleanUpAlg.cxx
4 /// \brief Remove raw digits corresponding to muon hits.
5 ///
6 /// \version $Id: TrackCleanUpAlg.cxx,v 1.11 2012-10-26 22:19:03 bckhouse Exp $
7 /// \author Kanika Sachdev (ksachdev@physics.umn.edu)
8 ////////////////////////////////////////////////////////////////////////
9 
11 
13 
14 #include "Geometry/Geometry.h"
16 #include "RecoBase/CellHit.h"
17 #include "RecoBase/RecoHit.h"
18 #include "RecoBase/WeightedHit.h"
19 #include "RecoBase/Cluster.h"
20 #include "RecoBase/Track.h"
21 #include "Utilities/func/MathUtil.h"
22 
23 #include "TGeoManager.h"
24 #include "TGeoNode.h"
25 
26 namespace murem
27 {
28 
29  /*********************************************************************************************/
30  ////////////////////////////////////////////////////////////////////////
31 
33  {
34  this->reconfigure(pset);
35  }
36 
37  /*********************************************************************************************/
38  ////////////////////////////////////////////////////////////////////////
39 
40 
42  {
43 
44  }
45 
46  /*********************************************************************************************/
47  ////////////////////////////////////////////////////////////////////////
48 
50  {
51  fMipRangeHigh = pset.get< float >("MipRangeHigh");
52  fMipRangeLow = pset.get< float >("MipRangeLow");
53  fMipValue = pset.get< float >("MipValue");
54  fVertexRegionDeDxCutOff = pset.get< float >("VertexRegionDeDxCutOff");
55 
56  fMipRangeHighInGev = pset.get< float >("MipRangeHighInGev");
57  fMipRangeLowInGev = pset.get< float >("MipRangeLowInGev");
58  fMipValueInGev = pset.get< float >("MipValueInGev");
59  fVertexRegionDeDxCutOffInGev = pset.get< float >("VertexRegionDeDxCutOffInGev");
60  }
61 
62  /*********************************************************************************************/
63  ////////////////////////////////////////////////////////////////////////
64 
65  /// Track Clean up close to the vertex.
66  /// For hits on the track, checks if the energy per length
67  /// in the first few planes is consistent with a MIP.
68  /// If it is not, modifies the hits and or the number of hits
69  /// associated with the track in a given plane.
70 
71  std::vector<rb::WeightedHit> TrackCleanUpAlg::CleanUpTrack(const rb::Track &muonTrack,
72  const rb::Cluster &slice,
73  double* mipRangeHigh,
74  double* mipRangeLow,
75  double* mipValue,
76  double* vertexRegionDeDxCutOff)
77  {
78 
79  // Check if the mip numbers are provided to us. If not, use
80  // defaults which are good only for muons.
81 
82  if(mipRangeHigh)
83  fMipRangeHigh = *mipRangeHigh;
84  if(mipRangeLow)
85  fMipRangeLow = *mipRangeLow;
86  if(mipValue)
87  fMipValue = *mipValue;
88 
89  // NB: everywhere in this script energy is calculated in units of MIP
90 
91  fExtraEOnTrack = 0;
92  fExtraEOnTrackPlanes.clear();
93 
94  double mipDeDx = fMipValue;
95  bool xViewCleaned = false;
96  bool yViewCleaned = false;
97  int vertexRegion[2];
98 
99  //if(vertexRegionDeDxCutOff)
100 
101  ComputeVertexRegion(muonTrack, vertexRegion, vertexRegionDeDxCutOff);
102 
103  int nPlanes = fPlanes.size();
104 
105  std::vector<rb::WeightedHit > trackHits;
106 
107  for(int i = 0; i < (int)muonTrack.NCell(); i++){
108  rb::WeightedHit newHit( muonTrack.Cell(i), 1 );
109  trackHits.push_back( newHit );
110  }
111 
112  for( int plidx = 0; plidx < nPlanes; plidx++ ){
113  if ( xViewCleaned && yViewCleaned )
114  break;
115 
116  int iPlane = fPlanes[plidx];
117  geo::View_t view = fgeom->Plane( iPlane )->View();
118 
119  if ( view == geo::kX){
120  if (iPlane > vertexRegion[0] ){
121  xViewCleaned = true;
122  continue;
123  }
124  }
125  else{
126  if ( iPlane > vertexRegion[1] ){
127  yViewCleaned = true;
128  continue;
129  }
130  }
131 
132  // Track clean is only done for the part of the track impacted
133  // by vertex energy. There can be variations in the dEdX beyond
134  // the MIP RMS due to delta rays escaping the cell etc. We want
135  // to make sure we don't have to deal with that by restricting
136  // ourselves to the first half.
137 
138  // This is the energy a MIP should deposit in this plane.
139  double eDep = (mipDeDx * fPlaneTrkLength[iPlane]);
140 
141  int flag = IsMIP( fPlaneDeDx[iPlane] );
142 
143  if(flag == 2){ //More than MIP
144 
145  SortByDistFromTrack( fPlaneHits[iPlane], muonTrack, iPlane );
146 
147  float trialE = 0;
148  int nPlaneHits = fPlaneHits[iPlane].size();
149  for( int iHit = 0; iHit< nPlaneHits; iHit++){
150  geo::OfflineChan p( iPlane, fPlaneHits[iPlane].at(iHit)->Cell() ) ;
151  trialE += fCellEnergy[p];
152 
153  int cellFlag = IsMIP( trialE/fPlaneTrkLength[iPlane]) ;
154 
155  if (cellFlag == 0){ // Not enough energy to make a MIP yet.
156 
157  rb::WeightedHit newHit( fPlaneHits[iPlane].at(iHit), 1);
158  fCleanUpWeightedHits.push_back(newHit);
159 
160  continue;
161  }
162  else{
163  // If we are not less than MIP, erase all
164  // other plane hits from trackHits
165 
166  for( int k = iHit+1; k<nPlaneHits; k++){
167  float weight = 0;
168  trackHits = ResetHitWeight( trackHits, fPlaneHits[iPlane].at(k), weight );
169  geo::OfflineChan pa( iPlane, fPlaneHits[iPlane].at(k)->Cell() ) ;
170  rb::WeightedHit newHit( fPlaneHits[iPlane].at(k), weight);
171  fCleanUpWeightedHits.push_back(newHit);
173  fExtraEOnTrackPlanes[iPlane] += fCellEnergy[pa];
174 
175  }
176 
177  if (cellFlag == 2){
178  double muEnergyInCell = eDep + fCellEnergy[p] - trialE;
179  float weight = muEnergyInCell/fCellEnergy[p];
180  trackHits = ResetHitWeight( trackHits, fPlaneHits[iPlane].at(iHit), weight );
181  rb::WeightedHit newHit(fPlaneHits[iPlane].at(iHit), weight);
182  fCleanUpWeightedHits.push_back(newHit);
183  fExtraEOnTrack += fCellEnergy[p] - muEnergyInCell;
184  fExtraEOnTrackPlanes[iPlane] += fCellEnergy[p] - muEnergyInCell;
185 
186  }
187 
188  break;
189  }
190 
191  }// end loop over fPlaneHits
192  }// end if on energy in plane is greater than MIP
193 
194 
195 
196  if( flag == 0){ //Less than MIP
197 
198  int nPlaneHits = fPlaneHits[iPlane].size();
199  unsigned int cell[ nPlaneHits ];
200  std::vector< art::Ptr<rb::CellHit> > slicePlaneHits;
201 
202  int nCells = slice.NCell();
203 
204  for( int h = 0; h < nPlaneHits; h++)
205  cell[h] = fPlaneHits[iPlane].at(h)->Cell();
206 
207  for( int hit = 0; hit < nCells; hit++){
208  for(int iHit = 0; iHit < nPlaneHits; iHit++){
209  if (slice.Cell(hit)->Plane() == iPlane &&
210  !(*(slice.Cell(hit)) == *( fPlaneHits[iPlane].at(iHit) )) )
211  slicePlaneHits.push_back(slice.Cell(hit));
212  }
213  }// end loop over slice hits
214 
215  SortByDistFromTrack( slicePlaneHits, muonTrack, iPlane );
216  float trialE = fPlaneEnergy[iPlane];
217 
218  int nSlicePlaneHits = slicePlaneHits.size();
219  for(int hit = 0; hit< nSlicePlaneHits; hit++){
220 
221  bool closeEnough = false;
222  for(int c = 0; c < nPlaneHits; c++ ){
223  int diff = slicePlaneHits[hit]->Cell() - cell[c];
224  if ( diff < 5 && diff > -5 )
225  closeEnough = true;
226  }
227 
228  if( closeEnough){
229 
230  // Need to add each slicePlaneHit to muonTrack so that it is
231  // calibrated accordingly
232 
233  rb::RecoHit rHit = muonTrack.RecoHit(slicePlaneHits[hit] );
234 
235  if( rHit.IsCalibrated() ){
236  trialE += rHit.MIP();
237 
238  int cellFlag = IsMIP( trialE/fPlaneTrkLength[iPlane] ) ;
239  geo::OfflineChan p( iPlane, slicePlaneHits[hit]->Cell() ) ;
240 
241  if (cellFlag == 0){
242  rb::WeightedHit newHit(slicePlaneHits[hit], 1);
243 
244  trackHits.push_back(newHit);
245  fCleanUpWeightedHits.push_back(newHit);
247  fExtraEOnTrackPlanes[iPlane] -= fCellEnergy[p];
248  }
249 
250  if (cellFlag == 1) {
251  rb::WeightedHit newHit(slicePlaneHits[hit], 1);
252 
253  trackHits.push_back(newHit);
254  fCleanUpWeightedHits.push_back(newHit);
256  fExtraEOnTrackPlanes[iPlane] -= fCellEnergy[p];
257  break;
258  }
259 
260  else if (cellFlag == 2){
261 
262  float tempEnergy = eDep + rHit.MIP() - trialE;
263 
264  rb::WeightedHit tempHit(slicePlaneHits[hit], tempEnergy/rHit.MIP() );
265  fExtraEOnTrack -= tempEnergy;
266  fExtraEOnTrackPlanes[iPlane] -= tempEnergy;
267  trackHits.push_back( tempHit );
268  break;
269  }
270 
271 
272  }
273  }// end of closeEnough
274  }// end loop over slice hits in iPlane
275  }// end if energy in plane is less than MIP
276 
277 
278  }// end loop over iPlanes.
279 
280  return trackHits;
281  }// end function CleanUpTrack
282 
283 
284 
285  /*********************************************************************************************/
286  ////////////////////////////////////////////////////////////////////////
287 
288  /// Checks if the given dE/dx is consistent with a MIP.
289  /// Returns 1 if dE/dx is consistent with MIP.
290  /// Returns 0, if dE/dx is less than MIP.
291  /// Returns 2, if dE/dx is more than MIP.
292 
293  int TrackCleanUpAlg::IsMIP( double dEdX ) const
294  {
295  // This will probably be changed to compare to MIP scale as
296  // returned by recohit, please put up with the hard coded
297  // numbers for now! Refer NOVA-doc-7658-v2, slide 11 for origin.
298  if ( dEdX < fMipRangeLow ) return 0;
299  // Less than MIP
300  if ( dEdX >= fMipRangeLow && dEdX <= fMipRangeHigh ) return 1;
301  // Consistent with MIP
302  if ( dEdX > fMipRangeHigh ) return 2;
303  // More than MIP
304  else return -1;
305  }//end function IsMIP
306 
307 
308  /// Finds and stores information plane by plane for a track.
309  /// The information it stores is energy per plane,
310  /// track length in plane, dE/dx in plane, energy per cell
311  /// and calibration constant for every cell on track.
312 
313 
314  std::map< unsigned int, double> TrackCleanUpAlg::DeDxInPlane(const rb::Track &muonTrack)
315  {
316 
317  fPlaneEnergy.clear();
318  fPlaneTrkLength.clear();
319  fPlaneDeDx.clear();
320  fPlanes.clear();
321  fPlanesX.clear();
322  fPlanesY.clear();
323  fPlaneHits.clear();
324 
325  int minPlane = muonTrack.MinPlane();
326  int maxPlane = muonTrack.MaxPlane();
327 
328  art::PtrVector<rb::CellHit> trkHits = muonTrack.AllCells();
329  SortByPlaneAndCell(trkHits);
330 
331 
332  // Check if the track is vertical. It's all very easy in that case
333  if(minPlane == maxPlane){
334  fPlanes.push_back(minPlane);
335  if( fgeom->Plane( minPlane)->View() == geo::kX)
336  fPlanesX.push_back( minPlane );
337  else
338  fPlanesY.push_back( maxPlane );
339 
340  double planeWidth = 2.0*fgeom->Plane(minPlane)->Cell(0)->HalfD();
341  // MIP() is calculated in one cell in z, calculate length in terms of cell width
342  fPlaneTrkLength[minPlane] = muonTrack.TotalLength()/planeWidth;
343  fPlaneEnergy[minPlane] = TrackEinMIP(muonTrack);
344  fPlaneDeDx[minPlane] = fPlaneEnergy[minPlane]/fPlaneTrkLength[minPlane];
345  for( auto const& iHit : trkHits){
346  fPlaneHits[ minPlane ].push_back( iHit );
347  rb::RecoHit rhit = muonTrack.RecoHit(iHit);
348  if( rhit.IsCalibrated() )
349  fCellEnergy[ iHit->OfflineChan()] = rhit.MIP();
350  }
351  }
352 
353  else{
354 
355 
356 
357  // Find all the unique planes that the track passes through
358  for( auto const& iHit : trkHits){
359  int iPlane = iHit->Plane();
360  fPlaneHits[ iPlane ].push_back( iHit );
361  rb::RecoHit rhit = muonTrack.RecoHit(iHit);
362  if( rhit.IsCalibrated() )
363  fCellEnergy[ iHit->OfflineChan()] = rhit.MIP();
364 
365  if( std::find( fPlanes.begin(), fPlanes.end(), iPlane) == fPlanes.end() )
366  // This plane is a new one!
367  fPlanes.push_back( iPlane );
368  if( fgeom->Plane( iPlane)->View() == geo::kX)
369  fPlanesX.push_back( iPlane );
370  else
371  fPlanesY.push_back( iPlane );
372 
373  if( rhit.IsCalibrated() )
374  fPlaneEnergy[ iPlane ] += rhit.MIP();
375  }// end loop over track hits
376 
377  // Now loop through the planes, find path-length through each and
378  // store the dE/dx.
379 
380 
381  for(int iPl = 0, nPlanes = fPlanes.size(); iPl < nPlanes; ++iPl){
382 
383  int iPlane = fPlanes[iPl];
384 
385  // We can give exaggerated results for vertical or close to vertical
386  // tracks, so just check the z direction cosine
387 
388  if( muonTrack.Dir().Z() < 0.2 && muonTrack.Dir().Z() > -0.2 ){
389  // Looks like we have a vertical-ish track on our hands.
390  // We will find the max and min cells in
391  // double cellWidth = 2*fgeom->Plane(iPlane)->Cell(0)->HalfD();
392  int ncell = fPlaneHits[ iPlane ].back()->Cell() - fPlaneHits[ iPlane ].front()->Cell();
393  // MIP() is calculated in one cell in z, calculate length in terms of cell width
394  fPlaneTrkLength[iPlane] = ncell - 1; // was cellWidth*( ncell - 1) previously
395  // the minus 1 here is assuming that on average only
396  // the half os the first and last cells in the plane
397  // are traversed.
398  }
399  else{
400  double planeHalfWidth = fgeom->Plane(iPlane)->Cell(0)->HalfD();
401  double cellCenter[3];
402  fgeom->Plane( iPlane )->Cell( 10 )->GetCenter(cellCenter);
403  double xbound1, ybound1, xbound2, ybound2;
404  muonTrack.InterpolateXY( cellCenter[2]-planeHalfWidth, xbound1, ybound1);
405  muonTrack.InterpolateXY( cellCenter[2]+planeHalfWidth, xbound2, ybound2);
406  TVector3 bound1( xbound1, ybound1, cellCenter[2]-planeHalfWidth);
407  TVector3 bound2( xbound2, ybound2, cellCenter[2]+planeHalfWidth);
408 
409  // MIP() is calculated in one cell in z, calculate length in terms of cell width
410  fPlaneTrkLength[iPlane] = (bound2 - bound1).Mag()/planeHalfWidth/2.0;
411  }
412 
413  fPlaneDeDx[iPlane] = fPlaneEnergy[iPlane]/fPlaneTrkLength[iPlane];
414  }// end loop over planes
415 
416  }
417 
418  return fPlaneDeDx;
419  }//End of DeDxInPlane
420 
421 
422  /*********************************************************************************************/
423  ////////////////////////////////////////////////////////////////////////
424 
425  /// Erases a specified hit from a std::vector of rb::CellHit.
426 
427  std::vector<rb::WeightedHit > TrackCleanUpAlg::ResetHitWeight( std::vector< rb::WeightedHit> &trackHits,
429  float &weight)
430  {
431  int nTrackHits = trackHits.size();
432  for(int iHit = 0; iHit < nTrackHits ; iHit++){
433  if(*(trackHits[iHit].hit) == *(hit))
434  trackHits[iHit].weight = weight;
435  }
436  return(trackHits);
437  }// End of EraseHit
438 
439 
440 
441  /*********************************************************************************************/
442  ////////////////////////////////////////////////////////////////////////
443 
444  // A function to compute the extent of hadronic energy impacted region
445  // along the length of the muon track.
446 
448  int *vertexRegion,
449  double* vertexRegionDeDxCutOff)
450  {
451  if(vertexRegionDeDxCutOff)
452  fVertexRegionDeDxCutOff = *vertexRegionDeDxCutOff;
453 
454  if( !(muonTrack.NCell() > 0) )
455  return;
456 
457  fAvgDeDx.clear();
458  DeDxInPlane(muonTrack);
459 
460  vertexRegion[0] = -1;
461  vertexRegion[1] = -1;
462 
463 
464  for( int view = 0; view < 2 ; view++){
465  std::vector<int> planes;
466 
467  if ( view == 0 )
468  planes = fPlanesX;
469  else
470  planes = fPlanesY;
471 
472  int nPlanes = planes.size();
473 
474  for(int ipl = 0; ipl < nPlanes - 2 ; ipl++ ){
475  float avgDeDx = 0;
476  for(int jpl = ipl; jpl < ipl+3; jpl++){
477  avgDeDx += fPlaneDeDx[ planes[jpl] ];
478  }
479 
480  fAvgDeDx[ planes[ipl] ] = avgDeDx/3;
481  }
482 
483  if( nPlanes > 7){
484  for(int ipl = 0; ipl < nPlanes - 3; ipl++ ){
485 
486  if( fAvgDeDx[ planes[ipl] + 1] < fVertexRegionDeDxCutOff &&
487  fAvgDeDx[ planes[ipl] + 2] < fVertexRegionDeDxCutOff &&
488  fAvgDeDx[ planes[ipl] + 3] < fVertexRegionDeDxCutOff &&
489  fAvgDeDx[ planes[ipl] + 4] < fVertexRegionDeDxCutOff){
490  vertexRegion[view] = planes[ipl];
491  break;
492  }
493  }
494  if( vertexRegion[view] == -1 && nPlanes>0 )
495  vertexRegion[view] = planes[ nPlanes - 1 ];
496  }
497  if( nPlanes < 8){
498  for(int ipl = 0; ipl < nPlanes - 3; ipl++ ){
499 
500  if( fPlaneDeDx[ planes[ipl] + 1] < fVertexRegionDeDxCutOff &&
501  fPlaneDeDx[ planes[ipl] + 2] < fVertexRegionDeDxCutOff &&
502  fPlaneDeDx[ planes[ipl] + 3] < fVertexRegionDeDxCutOff ){
503  vertexRegion[view] = planes[ipl];
504  break;
505  }
506  }
507  if( vertexRegion[view] == -1 && nPlanes>0 )
508  vertexRegion[view] = planes[ nPlanes - 1 ];
509  }
510  }
511  }// End of ComputeVertexRegion
512 
513 
514 
515  /*********************************************************************************************/
516  ////////////////////////////////////////////////////////////////////////
517 
518 
520  const rb::Track &track, int plane)
521  {
522 
523  int nHits = hits.size();
524 
525  for( int iHit = 0; iHit < nHits ; iHit++){
526 
527  TVector3 cellCenter, trkCenter;
528  fgeom->Plane( hits[iHit]->Plane() )->Cell( hits[iHit]->Cell() )->GetCenter( cellCenter );
529 
530  // Find the z position of this plane, and use rb::Track::InterpolateXY to
531  // figure the x,y position in this plane that the track passed through.
532  fgeom->Plane( plane )->Cell( hits[iHit]->Cell() )->GetCenter( trkCenter );
533  double x,y;
534  track.InterpolateXY( trkCenter[2], x, y );
535  TVector3 trkPoint( x, y, trkCenter[2]);
536 
537  float iCellDist = ( cellCenter - trkPoint ).Mag();
538 
539  for( int jHit = iHit+1; jHit < nHits; jHit++){
540  fgeom->Plane( hits[jHit]->Plane() )->Cell( hits[jHit]->Cell() )->GetCenter( cellCenter );
541 
542  float jCellDist = ( cellCenter - trkPoint ).Mag();
543 
544  if( jCellDist < iCellDist ){
545  art::Ptr<rb::CellHit> tmpCellHit = hits[iHit];
546  hits[iHit] = hits[jHit];
547  hits[jHit] = tmpCellHit;
548  iCellDist = jCellDist;
549  }
550  }//end inner loop to order hits in distance from track
551  }//end outer loop to order hits in distance from track
552  return;
553  }
554 
555 
556  /*********************************************************************************************/
557  ////////////////////////////////////////////////////////////////////////
558 
559 
561  const rb::Cluster &slice,
562  double* mipRangeHigh,
563  double* mipRangeLow,
564  double* mipValue,
565  double* vertexRegionDeDxCutOff)
566  {
567  CleanUpTrack( track, slice, mipRangeHigh, mipRangeLow, mipValue, vertexRegionDeDxCutOff);
568  return fExtraEOnTrack;
569  }
570 
572  const rb::Cluster &slice,
573  double* mipRangeHigh,
574  double* mipRangeLow,
575  double* mipValue,
576  double* vertexRegionDeDxCutOff)
577  {
578  std::cout<<"NB: entire TrackCleanUpAlg works in MIP units now, so your mipRange* values will be translated to the MIPs"<<std::endl;
579  std::cout<<"NB: if you haven't set mip arguments for this function default values from TrackCleanUpAlg.fcl will be called"<<std::endl;
580 
581  if(mipRangeHigh)
582  fMipRangeHighInGev = *mipRangeHigh;
583  if(mipRangeLow)
584  fMipRangeLowInGev = *mipRangeLow;
585  if(mipValue)
586  fMipValueInGev = *mipValue;
587  if(vertexRegionDeDxCutOff)
588  fVertexRegionDeDxCutOffInGev = *vertexRegionDeDxCutOff;
589 
590  double highInMIPs = fMipRangeHighInGev / fMipValueInGev;
591  double lowInMIPs = fMipRangeLowInGev / fMipValueInGev;
592  double mipInMIPs = 1.0;
593  double vertexRegionInMIPs = fVertexRegionDeDxCutOffInGev / fMipValueInGev;
594 
595  CleanUpTrack( track, slice, &highInMIPs, &lowInMIPs, &mipInMIPs, &vertexRegionInMIPs); // External packages will call this function with mipRange* in GeV/cm. Need to translate into MIP units
596  return fExtraEOnTrack/94.62; // This a constant in Calibrator_service.cc, GetGeVToMIPScale function
597  }
598 
599  /*********************************************************************************************/
600  ////////////////////////////////////////////////////////////////////////
601 
603  const rb::Cluster &slice,
604  double* mipRangeHigh,
605  double* mipRangeLow,
606  double* mipValue,
607  double* vertexRegionDeDxCutOff)
608  {
609 
610  CleanUpTrack( track, slice, mipRangeHigh, mipRangeLow, mipValue, vertexRegionDeDxCutOff);
611  return fExtraEOnTrackPlanes;
612  }
613 
615  const rb::Cluster &slice,
616  double* mipRangeHigh,
617  double* mipRangeLow,
618  double* mipValue,
619  double* vertexRegionDeDxCutOff)
620  {
621  std::cout<<"NB: entire TrackCleanUpAlg works in MIP units now, so your mipRange* values will be translated to the MIPs"<<std::endl;
622  std::cout<<"NB: if you haven't set mip arguments for this function default values from TrackCleanUpAlg.fcl will be called"<<std::endl;
623 
624  if(mipRangeHigh)
625  fMipRangeHighInGev = *mipRangeHigh;
626  if(mipRangeLow)
627  fMipRangeLowInGev = *mipRangeLow;
628  if(mipValue)
629  fMipValueInGev = *mipValue;
630  if(vertexRegionDeDxCutOff)
631  fVertexRegionDeDxCutOffInGev = *vertexRegionDeDxCutOff;
632 
633  double highInMIPs = fMipRangeHighInGev / fMipValueInGev;
634  double lowInMIPs = fMipRangeLowInGev / fMipValueInGev;
635  double mipInMIPs = 1.0;
636  double vertexRegionInMIPs = fVertexRegionDeDxCutOffInGev / fMipValueInGev;
637 
638  CleanUpTrack( track, slice, &highInMIPs, &lowInMIPs, &mipInMIPs, &vertexRegionInMIPs); // External packages will call this function with mipRange* in GeV/cm. Need to translate into MIP units
639  std::map<int, float> fExtraEOnTrackPlanesInGeV;
640  for (auto itr = fExtraEOnTrackPlanes.begin(); itr != fExtraEOnTrackPlanes.end(); itr++){
641  fExtraEOnTrackPlanesInGeV.insert ( std::pair<int, float>(itr->first, itr->second / 94.62) ); // This a constant in Calibrator_service.cc, GetGeVToMIPScale function
642  }
643  return fExtraEOnTrackPlanesInGeV;
644  }
645  /*********************************************************************************************/
646  ////////////////////////////////////////////////////////////////////////
647 
649  int *vertexPlanes,
650  double* vertexRegionDeDxCutOff)
651  {
652  int vertexRegion[2];
653  ComputeVertexRegion( track, vertexRegion, vertexRegionDeDxCutOff );
654  vertexPlanes[0] = vertexRegion[0];
655  vertexPlanes[1] = vertexRegion[1] ;
656  }
657 
658  /*********************************************************************************************/
659  ////////////////////////////////////////////////////////////////////////
660  std::vector<rb::WeightedHit> TrackCleanUpAlg::CleanUpWeightedHits(const rb::Track &track,
661  const rb::Cluster &slice,
662  double* mipRangeHigh,
663  double* mipRangeLow,
664  double* mipValue,
665  double* vertexRegionDeDxCutOff)
666  {
667 
668  CleanUpTrack( track, slice, mipRangeHigh, mipRangeLow, mipValue, vertexRegionDeDxCutOff);
669  return fCleanUpWeightedHits;
670  }
671 
672  /*********************************************************************************************/
673  ////////////////////////////////////////////////////////////////////////
675  {
676  double ret = 0;
677  for(unsigned int i = 0; i < track.NCell(); ++i){
678  const rb::RecoHit rhit = track.RecoHit(i);
679  if(rhit.IsCalibrated()) ret += rhit.MIP()*track.Weight(i);
680  }
681  return ret;
682 
683  }
684 
685 }// End of namespace murem
std::map< int, float > fExtraEOnTrackPlanes
A map of planes in the vertex region of a track and the extra energy found on them.
double HalfD() const
Definition: CellGeo.cxx:205
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
TrackCleanUpAlg(fhicl::ParameterSet const &pset)
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
tracking algorithms
std::map< int, std::vector< art::Ptr< rb::CellHit > > > fPlaneHits
A map to store the cell numbers that are in each plane of the track.
std::map< geo::OfflineChan, double > fCellEnergy
A map of cell and energy in Plane left by the track.
const Var weight
void reconfigure(const fhicl::ParameterSet &pset)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Definition: event.h:19
std::map< unsigned int, double > fPlaneDeDx
A map of plane and dE/dX of the track in the plane.
A collection of associated CellHits.
Definition: Cluster.h:47
std::vector< rb::WeightedHit > fCleanUpWeightedHits
A vector of rb::WeightedHit that are impacted by CleanUp procedure.
float fMipRangeHigh
Default upper bound of allowed MIP range (FHiCL parameter)
art::ServiceHandle< geo::Geometry > fgeom
Handle to Geometry service.
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
std::vector< rb::WeightedHit > ResetHitWeight(std::vector< rb::WeightedHit > &trackHits, art::Ptr< rb::CellHit > &hit, float &weight)
A function to reset the weight of a rb::WeightedHit in a given vector of WeightedHits.
const PlaneGeo * Plane(unsigned int i) const
float fVertexRegionDeDxCutOffInGev
the same as fVertexRegionDeDxCutOff but in GeV/cm units (FHiCL parameter)
void ComputeVertexRegion(const rb::Track &muonTrack, int *vertexRegion, double *vertexRegionDeDxCutOff=NULL)
A function to decide how large a vertex region should be considered for clean up. vertexRegion is an ...
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
std::map< int, float > ExtraEOnTrackPlanes(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
Calculates the total hadronic energy on the provided track based on the assumption of MIP for each pl...
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
int IsMIP(double dEdX) const
A function to check if the given dE/dx is consistent with MIP.
std::map< int, float > ExtraEOnTrackPlanesInGeV(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
float fVertexRegionDeDxCutOff
Default vertex region dedx cut off value (FHiCL parameter)
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
void hits()
Definition: readHits.C:15
float ExtraEOnTrack(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
Calculates the total hadronic energy on the provided track based on the assumption of MIP...
float fMipValueInGev
the same as fMipValue but in GeV/cm units (FHiCL parameter)
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
float fMipValue
Default MIP value in GeV/cm (FHiCL parameter)
float TrackEinMIP(const rb::Track &track)
A function for total energy calculation in terms of MIP.
std::map< unsigned int, double > DeDxInPlane(const rb::Track &muonTrack)
A function to calculate the energy, dE/dx and tracklength of a track through a plane.
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
std::vector< int > fPlanes
A vector of planes on muon track with hits in them.
OStream cout
Definition: OStream.cxx:6
std::vector< rb::WeightedHit > CleanUpTrack(const rb::Track &muonTrack, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
A function to disentangle the muon track from the hadronic energy close to the vertex. Returns a vector of rb::WeightedHits where the hit is the cellhits that were on the given track and weight is the fraction of the energy in the hit that came from the track. mipRangeLow, mipRangeHigh, mipValue and vertexRegionDeDxCutOff are pointers that point to the value of mip range limits, peak values andthe dedx cut off used to find the vertex region. If these pointers are not provided, the default values are used. The default values are good for muons, not for other particles.
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
float fExtraEOnTrack
FHiCL Parameters.
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
float MIP() const
Definition: RecoHit.cxx:58
std::vector< rb::WeightedHit > CleanUpWeightedHits(const rb::Track &muonTrack, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
Returns a vector of rb::WeightedHit that have been impacted by the CleanUp procedure. hits are all the hits on the muonTrack within the vertex region and weight indicates the fraction of energy deposited in the cell that belonged to the track. mipRangeLow, mipRangeHigh, mipValue and vertexRegionDeDxCutOff are pointers that point to the value of mip range limits, peak values and the dedx cut off used to find the vertex region. If these pointers are not provided, the default values are used. The default values are good for muons, not for other particles.
std::vector< int > fPlanesX
A vector of X-view planes on muon track with hits in them.
Definition: event.h:1
A (plane, cell) pair.
Definition: OfflineChan.h:17
std::vector< int > fPlanesY
A vector of Y-view planes on muon track with hits in them.
float Mag() const
Simple hit+weight pair, returned from rb::Cluster::WeightedHits.
Definition: WeightedHit.h:18
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
virtual void InterpolateXY(double z, double &x, double &y) const
Definition: Track.cxx:325
void SortByDistFromTrack(std::vector< art::Ptr< rb::CellHit > > &hits, const rb::Track &track, int plane)
Sorts a vector of cellhits in ascending order by distance from the trajPt on track.
void SortByPlaneAndCell(std::vector< rb::WeightedHit > &whits)
std::map< unsigned int, double > fPlaneTrkLength
A map of plane and track length through plane.
float fMipRangeLowInGev
the same as fMipRangeLow but in GeV/cm units (FHiCL parameter)
void LastVertexPlanes(const rb::Track &muonTrack, int *vertexPlanes, double *vertexRegionDeDxCutOff=NULL)
A function to decide how large a vertex region should be considered for clean up. vertexRegion is an ...
float ExtraEOnTrackInGeV(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
double Weight(unsigned int globalIdx) const
Weight assigned to the cell.
Definition: Cluster.cxx:209
Encapsulate the geometry of one entire detector (near, far, ndos)
std::map< unsigned int, double > fPlaneEnergy
A map of plane and energy in plane left by the track.
float fMipRangeHighInGev
the same as fMipRangeHigh but in GeV/cm units (FHiCL parameter)
std::map< unsigned int, double > fAvgDeDx
A map of plane and average dedx over this and the next 2 planes.
float fMipRangeLow
Default lower bound of allowed MIP range (FHiCL parameter)