DeconvolveAlg.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file DeconvolveAlg.cxx
3 /// \brief Separates the hadronic energy and electromagnetic energy in
4 /// a cell
5 ///
6 /// \author Kanika Sachdev (ksachdev@physics.umn.edu)
7 ////////////////////////////////////////////////////////////////////////
8 
12 
14 
15 #include "Calibrator/Calibrator.h"
17 #include "GeometryObjects/Geo.h"
18 #include "Utilities/func/MathUtil.h"
19 
20 #include "TGeoManager.h"
21 #include "TGeoNode.h"
22 #include "TMatrixD.h"
23 
24 #include <string.h>
25 
26 namespace slid
27 {
28 
29  /*********************************************************************************************/
30  ////////////////////////////////////////////////////////////////////////
31 
33  {
34  this->reconfigure(pset);
35  mf::LogInfo("DeconvolveAlg") << " DeconvolveAlg::DeconvolveAlg()\n";
36  }
37 
38  /*********************************************************************************************/
39  ////////////////////////////////////////////////////////////////////////
40 
42  {
43 
44  }
45 
46  /*********************************************************************************************/
47  ////////////////////////////////////////////////////////////////////////
48 
50  {
51  fExpConst = pset.get< float >("ExpConst");
52  }
53 
54  /*********************************************************************************************/
55  ////////////////////////////////////////////////////////////////////////
56 
57  void DeconvolveAlg::GetWeights( std::vector<rb::Prong> & prongs,
58  const std::vector< float > & totpe)
59  {
60  std::vector<rb::WeightedHit> wHits;
61  const int nPro = prongs.size();
62 
63  // This map holds info about the cells that are shared by
64  // multiple prongs. The key is the channel, and the value
65  // is another map. In the inner map, the key is the index
66  // of the prong and the value is the index of the cell in
67  // that prong. Therefore, proIdx[chan][prgIdx] = [cellIdx]
68  std::map< uint32_t , std::map< int, int > > proIdx ;
69 
70  for( int iPro = 0; iPro < nPro; ++iPro){
71 
72  // loop over cells in ith prong
73  for( unsigned int iCell = 0; iCell < prongs[iPro].NCell(); ++iCell ){
74  const art::Ptr<rb::CellHit>& ihit = prongs[iPro].Cell(iCell);
75 
76  // In the following loop, we will find which other
77  // prongs the ihits[iCell] belongs to.
78 
79  for( int jPro = iPro+1; jPro < nPro; ++jPro){
80 
81  // loop over cells in jth prong
82  for( unsigned int jCell = 0 ; jCell < prongs[jPro].NCell(); ++jCell){
83  const art::Ptr<rb::CellHit>& jhit = prongs[jPro].Cell(jCell);
84 
85  if( *jhit == *ihit){
86  proIdx[ ihit->Channel() ][ iPro ] = iCell;
87  proIdx[ ihit->Channel() ][ jPro ] = jCell;
88  }
89 
90  }// end loop over jprong cells
91  }// end j loop over prongs
92  }// end loop over iprong cells
93  }// end i loop over prongs
94 
95  // We now have all the info we need to set weights:
96  for( auto & cmap : proIdx){
97  double weights[ cmap.second.size()];
98  double total = 0;
99  int count = 0;
100  for( auto & pmap : cmap.second){
101 
102  rb::Prong iProng = prongs[ pmap.first ];
103  int iCell = pmap.second;
104  double dist = DistanceToCore( iProng.Cell( iCell), iProng );
105  if (dist == 9999)
106  weights[count] = -1;
107  else{
108 
109  // Shower lateral profile constant from a fit to
110  // (Transverse dE/dx)/(shower Energy) as a function
111  // of distance to shower core for Monte Carlo electron.
112 
113  weights[count] = totpe[ pmap.first] * exp( -fExpConst*dist);
114  total += weights[count];
115  }
116  count = count +1;
117  }// end loop over proIdx.second
118 
119  count=0;
120  for( auto & pmap : cmap.second){
121 
122  if( weights[count] >= 0 )
123  prongs[pmap.first].SetWeight( pmap.second, weights[count]/total);
124  else
125  prongs[pmap.first].SetWeight( pmap.second, 0.0 );
126  count = count +1;
127  }
128  }// end loop over proIdx
129 
130  return;
131 
132  }// end of GetWeights
133 
134 
135  /*********************************************************************************************/
136  ////////////////////////////////////////////////////////////////////////
137 
139  rb::Prong & prong)
140  {
141 
142  double dist = 9999, di;
143 
144  // If the prong is 2D, reclustering a cell from the other view
145  // makes no sense. Bail.
146  if( prong.Is2D() && ( prong.View() != cHit->View() ) )
147  return dist;
148 
149  TVector3 start = prong.Start();
150  TVector3 dir = prong.Dir();
151 
152  TVector3 xyz;
153  fGeom->Plane(cHit->Plane())->Cell(cHit->Cell())->GetCenter(xyz);
154  const TVector3 dxyz = (cHit->View() == geo::kX) ? TVector3(0, 1, 0): TVector3(1, 0, 0);
155 
156  // geo::Closest approach takes in 4 3D points-
157  // first two are start and stop points of a
158  // line and the nxt two, of the other.
159  // It returns the square of the distance of
160  // closest approach between two line.
161 
162  dist = geo::ClosestApproach(start, start+dir,
163  xyz, xyz+dxyz, &di);
164 
165  return sqrt(dist);
166 
167  } // end of DistanceToCore
168 
169 
170  /*********************************************************************************************/
171  ////////////////////////////////////////////////////////////////////////
172 
173 
174 
175  /*********************************************************************************************/
176  ////////////////////////////////////////////////////////////////////////
177 
178 }// End of namespace jmprong
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
void reconfigure(const fhicl::ParameterSet &pset)
Preforms energy deconvolution for cells that may belong to multple clusters. It does so...
bool Is2D() const
Definition: Cluster.h:96
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
float DistanceToCore(art::Ptr< rb::CellHit > cHit, rb::Prong &prong)
unsigned short Plane() const
Definition: CellHit.h:39
art::ServiceHandle< geo::Geometry > fGeom
Definition: DeconvolveAlg.h:54
DeconvolveAlg(fhicl::ParameterSet const &pset)
geo::View_t View() const
Definition: CellHit.h:41
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
double dist
Definition: runWimpSim.h:113
unsigned short Cell() const
Definition: CellHit.h:40
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Collect Geo headers and supply basic geometry functions.
float fExpConst
The constant length scale used in energy deconvolution.
Definition: DeconvolveAlg.h:55
Var weights
void GetWeights(std::vector< rb::Prong > &prong, const std::vector< float > &totpe)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
A Cluster with defined start position and direction.
Definition: Prong.h:19
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
Definition: Geo.cxx:222
TDirectory * dir
Definition: macro.C:5
uint32_t Channel() const
Definition: RawDigit.h:84
cmap::CMap class source code
Definition: CMap.cxx:17
Build slid::LID objects to store electron ID, if asked for, otherwise, calculate LID info and make av...
Definition: FillPIDs.h:13