SignalBasis.h 9.51 KB
Newer Older
1 2 3 4 5 6 7
#ifndef SIGNALBASIS_H_INCLUDED
#define SIGNALBASIS_H_INCLUDED

#include <memory.h>
#include <cmath>
#include <string>

dino's avatar
dino committed
8 9 10
#include "GridSearchParams.h"
#include "GridSearchClasses.h"

11 12 13 14 15 16 17 18
// As default, the fine grid is built using a cube surrounding the selected coarse point
// Alternatively it is possible to select one of the following
//#define ASSIGN_SPHERE         // multiply assigned points filling a surrounding sphere
//#define ASSIGN_TO_NEAREST     // points of second pass are assigned to only one coarse (not working well)

// Parameters of the coarse-fine grid
// Normally we work with the 2mm bases. The 1mm bases are too big and slow
// and are used occasionally just to check the behaviour of the 2mm ones.
19 20
//#define GRID1MM
#ifdef GRID1MM
21 22 23 24 25 26
  const int fstep  = 1;         // spacing (mm) of original grid
  const int cstep  = 3;         // spacing (mm) of coarse   grid
#if defined(ASSIGN_SPHERE) || defined(ASSIGN_NEAREST)
  const int cstepx = cstep;     // spacing (mm) of coarse   grid
  const int cstepy = cstep;     // spacing (mm) of coarse   grid
  const int cstepz = cstep;     // spacing (mm) of coarse   grid
27
#else
28 29 30 31 32 33
  const int cstepx = 3;         // spacing (mm) of coarse   grid
  const int cstepy = 3;         // spacing (mm) of coarse   grid
  const int cstepz = 3;         // spacing (mm) of coarse   grid
#endif
#else
  const int fstep  = 2;         // spacing (mm) of original grid
34
  const int cstep  = 6;         // spacing (mm) of coarse   grid
35 36 37 38 39
#if defined(ASSIGN_SPHERE) || defined(ASSIGN_NEAREST)
  const int cstepx = cstep;     // spacing (mm) of coarse   grid
  const int cstepy = cstep;     // spacing (mm) of coarse   grid
  const int cstepz = cstep;     // spacing (mm) of coarse   grid
#else
40 41 42
  const int cstepx = cstep;     // spacing (mm) of coarse   grid
  const int cstepy = cstep;     // spacing (mm) of coarse   grid
  const int cstepz = cstep;     // spacing (mm) of coarse   grid
43
#endif
44 45
#endif

46

47 48 49 50 51
// If the number of grid points in a segment is N, the optimum step of the coarse grid c is obtained
// by minimising M ~ N/c^3 + k*c^3  (coarse points + fine-points around it, k~4(sphere) k~8(cube))
// This gives c ~ pow(N/k, 1/6)
// With fstep=2,  N~1000 --> 2.5(k=4) ~2.2(k=8) ) with gains of ~5   ==> cstep = 2 (or 3)
// With fstep=1,  N~8000 --> 3.5(k=4) ~3.2(k=8) ) with gains of ~20  ==> cstep = 3
52

53
// Even if the calculated signals are almost always given in time steps of 5 ns,
54
// the basis in memory is always in "natural" steps of 10 ns.
55
// Even if the calculated signals are, by definition, zero before "T=0", we keep the code
56 57
// consistent with a small pretrigger (BZERO, usually set at 0). This is to allow sliding
// T0(calc) against the time-zero of the experimental traces (DZERO, normally set at 8).
58

59 60 61 62
class SignalBasis  
{
public:
  std::string theFileName;
63 64 65
  bool        bCenters;       // the center of segments has been calculated
  bool        bFullGrid;      // full-grid is in memory 
  bool        bCoarseFine;    // the coarse-fine structures are valid
66
  bool        bCoarseOnly;    // only the coarse part of the coarse-fine structures is valid
67
  int         numPts[NSEGS];  // number of points in the segments
dino's avatar
dino committed
68
  pointPsa   *segPts[NSEGS];  // the points of the segment, generated in Read()
69
  pointPsa    averPt[NCHAN];  // the "center" of the segment (in INDCC the overall center)
dino's avatar
dino committed
70
  float      *ddPtPt[NSEGS];  // precalculated point-point distance inside a segment
71 72 73 74

  float xTalkProp[NSEGS][NSEGS];
  float xTalkDiff[NSEGS][NSEGS];
  float xTalkAver[NSEGS];
75
  bool  segBroken[NSEGS+2];   // segBroken[NSEGS] always false; segBroken[NSEGS+1] true if any broken segment is present
76

77
  bool  xTalkCorr;
78

79 80 81 82 83 84
  // Bart's file format
  typedef struct {
    float T0;                 // time of interaction (in us)
    float Eint;               // energy fractions == 1
    float Pos[3];             // x y z (mm)
    int   Trig;               // Core trigger time (sample index)
85 86
    float Es[ 1+NSEGS     ];  // energy in segments (0=core)
    float Tr[(1+NSEGS)*120];  // CC+nSeg traces given in steps of 5 ns. The first 10 samples are empty (pretrigger of 50 ns)
87
  } bbhit_t;
88

89 90 91 92 93 94 95 96

  // Jeremie's file format
  typedef struct {
      int   SegTrig;              // Segment containing the interaction
      float Pos[3];		  // x y z (mm)
      float Tr[(NSEGS+1)*BSIZE];  // CC+nSeg traces given in steps of 10 ns. 
  } Newbbhit_t;

97
  // coarse grid
98 99 100 101 102 103 104 105
  typedef struct {
    float  x;   // coordinates the point
    float  y;
    float  z;
    float  d;
    int    j;   // which basis point it corresponds
    int    n;   // how many points of the original grid get associated with it
  } cgrid;
106

107
  // coarse-fine grid
108 109 110 111 112 113
  typedef struct {
    int  numCoarse;     // how many coarse grid points
    int  numFine;       // length of fine grid points list
    int *posCoarse;     // first of the fine grid indices in indFine for this coarse
    int *indFine;       // list of fine points
  } cfgrid;
114

115 116 117
  cfgrid cflist[NSEGS];

  SignalBasis() {
118
    bCenters = bFullGrid = bCoarseFine = bCoarseOnly = xTalkCorr = false;
119
    memset(numPts,    0,     sizeof(numPts));
dino's avatar
dino committed
120 121
    memset(segPts,    0,     sizeof(segPts));
    memset(ddPtPt,    0,     sizeof(ddPtPt));
122
    memset(cflist,    0,     sizeof(cflist));
123
    memset(segBroken, false, sizeof(segBroken));
124 125 126
  }
  virtual ~SignalBasis() {
    for(int ns=0; ns<NSEGS; ns++) {
dino's avatar
dino committed
127 128
      if(segPts[ns]) {
        delete[] segPts[ns]; segPts[ns] = NULL;
129
      }
130
      if(bCoarseFine) {
131 132
        delete [] cflist[ns].indFine;   cflist[ns].indFine   = NULL;
        delete [] cflist[ns].posCoarse; cflist[ns].posCoarse = NULL;
133 134 135
      }
    }
  }
dino's avatar
dino committed
136 137 138 139 140 141 142 143 144 145 146 147 148 149
  int   ReadBasis(std::string fname,    bool keep = true);
  int   ReadBasisFormatVentu(FILE * fp, bool keep);
  int   ReadBasisFormatBartB(FILE * fp, bool keep, int numSignals);
  int   ReadSegmentCenters(std::string fname);
  void  CalcPtPtDistance(int ns);
  float PtPtDistance(int ns, int n1, int n2);
  int   MakeCoarseFineList(bool both, bool verbose);
  bool  CoarseFineForOneSegment(int ns, bool both, bool verbose);
  bool  AssignIt(float px, float py, float pz, float cx, float cy, float cz);
  bool  IsCoarse(int jj, cgrid *cgr, int numC);
  int   NearestPoint(float px, float py, float pz, int ns, float& dist);
  int   NearestCoarse(float px, float py, float pz, cgrid *cgr, int numC, float& dist);
  void  SmoothBasisSignals(int nsmoo);
  void  ExperimentalResponse();
150
  void  ExponentialResponse(const float *TauSGCC, const float TauDecay, float smooth);
151
  void  CalculateSegmentCenters(bool verb);
dino's avatar
dino committed
152 153 154 155
  void  FindNeighbours();
  int   StoreNetCharges();
  void  CalculateSumOfSignals();
  int   DifferentiateNetChargeSignals(int lag);
156
  int   DifferentiateAllSignals(int lag);
dino's avatar
dino committed
157 158 159 160
  void  ReorderSamples();
  bool  ReadXtalkCoeffs(std::string cname, float xfactor, bool verbose);
  bool  CalcXtalkCoeffs(float *xTprop, float xfactor);
  bool  XtalkCorrection(int badSegment);
161 162
  float GetStepCoarse() {return float(cstep);}
  float GetStepFine()   {return float(fstep);}
163 164
};

dino's avatar
dino committed
165 166 167 168 169
inline float SignalBasis::PtPtDistance(int ns, int n1, int n2)
{
  return ddPtPt[ns][n1*numPts[ns] + n2];
}

170 171 172 173 174 175 176 177 178 179 180 181 182
inline bool SignalBasis::AssignIt(float px, float py, float pz, float cx, float cy, float cz)
{
  //const float surface = 0;              // to exclude the surface
  const float surface = 0.1f;           // to include the surface
  //const float surface = fstep + 0.1f;   // to extend the surface to the next layer
#ifdef ASSIGN_SPHERE
  // all points inside a sphere centered on this coarse grid point
  float dx = px - cx;
  float dy = py - cy;
  float dz = pz - cz;
  float d2 = dx*dx + dy*dy + dz*dz;
  bool inside = ( d2 < (cstep+surface)*(cstep+surface) );
#else
183
  // all points inside a "cube" centered on this coarse grid point
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201
  float dx = px - cx; if(dx < 0) dx = -dx;
  float dy = py - cy; if(dy < 0) dy = -dy;
  float dz = pz - cz; if(dz < 0) dz = -dz;
  bool inside = ( dx<(cstepx+surface) && dy<(cstepy+surface) && dz<(cstepz+surface) );
#endif
  return inside;
}

inline bool SignalBasis::IsCoarse(int jj, cgrid *cgr, int numC)
{
  for(int cc = 0; cc < numC; cc++) {
    if(cgr[cc].j == jj) {
      return true;
    }
  }
  return false;
}

202
// input value of dist used to stop search when a point closer than this
203 204
inline int SignalBasis::NearestPoint(float px, float py, float pz, int ns, float& dist)
{
205 206
  int       npts = numPts[ns];  // number of points in this segment
  pointPsa *ppts = segPts[ns];  // the points of this segment
207 208

  int   cnear = -1;
209 210
  float cprec = dist*dist;

211 212 213 214 215 216 217 218 219
  float cdist = (float)1.e20;
  for(int jj = 0; jj < npts; jj++) {
    float dx = ppts[jj].x - px;
    float dy = ppts[jj].y - py;
    float dz = ppts[jj].z - pz;
    float d2 = dx*dx + dy*dy + dz*dz;
    if(d2 < cdist) {
      cdist = d2;
      cnear = jj;
220 221
      if(cdist < cprec)
        break;
222 223 224 225 226 227
    }
  }
  dist = (float)sqrt(cdist);
  return cnear;
}

228
// input value of dist used to stop search when a point closer than this
229 230 231
inline int SignalBasis::NearestCoarse(float px, float py, float pz, cgrid *cgr, int numC, float& dist)
{
  int   cnear = -1;
232 233
  float cprec = dist*dist;

234 235 236 237 238 239 240 241 242
  float cdist = (float)1.e20;
  for(int kk = 0; kk < numC; kk++) {
    float dx = cgr[kk].x - px;
    float dy = cgr[kk].y - py;
    float dz = cgr[kk].z - pz;
    float d2 = dx*dx + dy*dy + dz*dz;
    if(d2 < cdist) {
      cdist = d2;
      cnear = kk;
243 244
      if(cdist < cprec)
        break;
245 246 247 248 249 250
    }
  }
  dist = (float)sqrt(cdist);
  return cnear;
}

251 252
#endif  // SIGNALBASIS_H_INCLUDED