DistanceMap.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file DistanceMap.cxx
3 /// \brief Calculate and cache electrostatic potential between cells
4 /// \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6
7 #include "LEM/func/DistanceMap.h"
8
9 #include "TH2.h"
10
11 #include <cstdlib>
12
13 #include "Utilities/func/MathUtil.h"
14
15
16 namespace lem
17 {
18  // Best aspect ratio empirically. Matches fcl default
19  double DistanceMap::fgPlaneScale = 3.5;
20  double DistanceMap::fgCellScale = 10.5;
21  // 1/sqrt(dr). Best power empirically. EM would be -1
22  double DistanceMap::fgDecayPower = -0.25;
23  bool DistanceMap::fgExpMode = false;
24
25  //......................................................................
27  {
28  static DistanceMap instance;
29  return instance;
30  }
31
32  //......................................................................
34  {
35  for(int p = 0; p < kNumPlanes; ++p){
36  for(int c = 0; c < kNumCells; ++c){
37  // For long enough distances don't bother to do the integral over cells
38  if(p > 5 || c > 5){
39  fInvDistMap[p][c] = DistFunc(p, c);
40  continue;
41  }
42
43  double E = 0;
44
45  const int kNumSteps = 10;
46  for(int dp1 = 0; dp1 < kNumSteps; ++dp1){
47  for(int dp2 = 0; dp2 < kNumSteps; ++dp2){
48  const double dp = (dp1-dp2)/double(kNumSteps);
49  for(int dc1 = 0; dc1 < kNumSteps; ++dc1){
50  for(int dc2 = 0; dc2 < kNumSteps; ++dc2){
51  // Avoid division by zero
52  if(p == 0 && c == 0 && dp1 == dp2 && dc1 == dc2) continue;
53
54  const double dc = (dc1-dc2)/double(kNumSteps);
55
56  E += DistFunc(p+dp, c+dc);
57  }
58  }
59  }
60  } // end for dp1
61
62  // Don't count the terms that we continue'd out of
63  int count = util::sqr(util::sqr(kNumSteps));
64  if(p == 0 && c == 0) count -= util::sqr(kNumSteps);
65
66  fInvDistMap[p][c] = E/count;
67  } // end for c
68  } // end for p
69  }
70
71  //......................................................................
72  TH2* DistanceMap::MakeHist() const
73  {
74  TH2* hDistMap = new TH2F("", ";#Deltaplanes;#Deltacells", 5, 0, 5, 5, 0, 5);
75
76  for(int p = 0; p < kNumPlanes; ++p){
77  for(int c = 0; c < kNumCells; ++c){
78  hDistMap->Fill(p, c, 1/fInvDistMap[p][c]);
79  } // end for c
80  } // end for p
81
82  return hDistMap;
83  }
84
85  //......................................................................
86  double DistanceMap::DistFunc(double dp, double dc)
87  {
88  const double dr = util::pythag(fgPlaneScale*dp, fgCellScale*dc);
89
90  if(fgExpMode){
91  // Can control the coefficient with planeScale and cellScale
92  return exp(-dr);
93  }
94
95  return pow(dr, fgDecayPower);
96  }
97
98  //......................................................................
99  double DistanceMap::InvDist(int aplane, int bplane, int acell, int bcell) const
100  {
101  const int dp = abs(aplane-bplane);
102  const int dc = abs(acell-bcell);
103
104  // Very far away
105  if(dp < 0) return 0;
106  if(dp >= kNumPlanes) return 0;
107  if(dc < 0) return 0;
108  if(dc >= kNumCells) return 0;
109
110  // Comparison between different views shouldn't happen
111  assert(dp%2 == 0);
112
113  return fInvDistMap[dp][dc];
114  }
115
116  //......................................................................
117  double DistanceMap::InvDist(const lem::LiteHit& a, const lem::LiteHit& b) const
118  {
119  return InvDist(a.Plane(), b.Plane(), a.Cell(), b.Cell());
120  }
121
122 } // namespace
static double DistFunc(double dp, double dc)
Definition: DistanceMap.cxx:86
const char * p
Definition: xmltok.h:285
constexpr T pow(T x)
Definition: pow.h:75
void abs(TH1 *hist)
double fInvDistMap[kNumPlanes][kNumCells]
Definition: DistanceMap.h:45
int Cell() const
Definition: LiteHit.h:39
static double fgDecayPower
Definition: DistanceMap.h:42
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
PID
Definition: FillPIDs.h:14
int Plane() const
Definition: LiteHit.h:38
static const DistanceMap & Instance()
Singleton.
Definition: DistanceMap.cxx:26
Float_t E
Definition: plot.C:20
const double a
TH2 * MakeHist() const
Definition: DistanceMap.cxx:72
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
static double fgCellScale
Definition: DistanceMap.h:41
static bool fgExpMode
Definition: DistanceMap.h:43
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
Compressed hit info, basic component of LEM events.
Definition: LiteHit.h:18
const int kNumCells
Definition: DistanceMap.h:17
Calculate and cache electrostatic potential between cells.
double InvDist(int aplane, int bplane, int acell, int bcell) const
Definition: DistanceMap.cxx:99
const int kNumPlanes
Definition: DistanceMap.h:16
static double fgPlaneScale
Definition: DistanceMap.h:40
Calculate and cache electrostatic potential between cells.
Definition: DistanceMap.h:20