Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DMaterialMap.cc
Go to the documentation of this file.
1 
2 
3 #include <iostream>
4 #include <set>
5 #include <cmath>
6 using namespace std;
7 
8 #include <DANA/DApplication.h>
9 #include <DVector2.h>
10 
11 #include "DMaterialMap.h"
12 #include "DMagneticFieldMap.h"
13 
14 
15 //-----------------
16 // DMaterialMap (Constructor)
17 //-----------------
18 DMaterialMap::DMaterialMap(string namepath, JCalibration *jcalib)
19 {
20  /// Read the specified material map in from the calibration database.
21  /// This will read in the map and figure out the number of grid
22  /// points in each direction (r, and z) and the range in each.
23 
24  IS_VALID = false;
25  MAX_BOUNDARY_SEARCH_STEPS = 30;
26  ENABLE_BOUNDARY_CHECK = true;
27 
28  gPARMS->SetDefaultParameter("GEOM:MAX_BOUNDARY_SEARCH_STEPS", MAX_BOUNDARY_SEARCH_STEPS, "Maximum number of steps (cells) to iterate when searching for a material boundary in DMaterialMap::EstimatedDistanceToBoundary(...)");
29  gPARMS->SetDefaultParameter("GEOM:ENABLE_BOUNDARY_CHECK", ENABLE_BOUNDARY_CHECK, "Enable boundary checking (superceeds any setting in DReferenceTrajectory). This is for debugging only.");
30 
31  this->namepath = namepath;
32 
33  // Read in map from calibration database. This should really be
34  // built into a run-dependent geometry framework, but for now
35  // we do it this way.
36  this->jcalib = jcalib;
37  if(!jcalib)return;
38 
39  string blank(' ',80);
40  //cout<<blank<<"\r"; cout.flush(); // clear line
41  //cout<<"Reading "<<namepath<<" ... "; cout.flush();
42  vector< vector<float> > Mmap;
43  jcalib->Get(namepath, Mmap);
44  //cout<<(int)Mmap.size()<<" entries (";
45  if(Mmap.size()<1){
46  //cout<<")"<<endl;
47  return;
48  }
49 
50  // The map should be on a grid with equal spacing in r, and z.
51  // Here we want to determine the number of points in each of these
52  // dimensions and the range.
53  // The easiest way to do this is to use a map<float, int> to make a
54  // histogram of the entries by using the key to hold the extent
55  // so that the number of entries will be equal to the number of
56  // different values.
57  map<float, int> rvals;
58  map<float, int> zvals;
59  double rmin, zmin, rmax, zmax;
60  rmin = zmin = 1.0E6;
61  rmax = zmax = -1.0E6;
62  for(unsigned int i=0; i<Mmap.size(); i++){
63  vector<float> &a = Mmap[i];
64  float &r = a[0];
65  float &z = a[1];
66 
67  rvals[r] = 1;
68  zvals[z] = 1;
69  if(r<rmin)rmin=r;
70  if(z<zmin)zmin=z;
71  if(r>rmax)rmax=r;
72  if(z>zmax)zmax=z;
73  }
74  Nr = rvals.size();
75  Nz = zvals.size();
76  r0 = rmin;
77  z0 = zmin;
78  dr = (rmax-rmin)/(double)(Nr-1);
79  dz = (zmax-zmin)/(double)(Nz-1);
80 
81  one_over_dr=1./dr;
82  one_over_dz=1./dz;
83  //cout<<"Nr="<<Nr;
84  //cout<<" Nz="<<Nz;
85  //cout<<")"<<endl;
86 
87  // The values in the table are stored with r,z at the center of the node.
88  // This means the actual map extends half a bin further out than the current
89  // (local variable) rmin,rmax and zmin,zmax values. Set the class data
90  // members to be the actual map limits
91  this->rmin = rmin-dr/2.0;
92  this->rmax = rmax+dr/2.0;
93  this->zmin = zmin-dz/2.0;
94  this->zmax = zmax+dz/2.0;
95 
96  // Set sizes of nested vectors to hold node data
97  nodes.resize(Nr);
98  for(int ir=0; ir<Nr; ir++){
99  nodes[ir].resize(Nz);
100  }
101 
102  // Fill table
103  for(unsigned int i=0; i<Mmap.size(); i++){
104  vector<float> &a = Mmap[i];
105  float &r = a[0];
106  float &z = a[1];
107  int ir = (int)floor((r-this->rmin)/dr);
108  int iz = (int)floor((z-this->zmin)/dz);
109  if(ir<0 || ir>=Nr){_DBG_<<"ir out of range: ir="<<ir<<" Nr="<<Nr<<endl; continue;}
110  if(iz<0 || iz>=Nz){_DBG_<<"iz out of range: iz="<<iz<<" Nz="<<Nz<<endl; continue;}
111  MaterialNode &node = nodes[ir][iz];
112  node.A = a[2];
113  node.Z = a[3];
114  node.Density = a[4];
115  node.RadLen = a[5];
116  node.rhoZ_overA = a[6];
117  node.rhoZ_overA_logI = a[7];
118  node.LogI=node.rhoZ_overA_logI/(node.rhoZ_overA==0.0 ? 1.0:node.rhoZ_overA);
119  node.KrhoZ_overA=0.1535e-3*node.rhoZ_overA;
120  node.chi2c_factor=a[8];
121  node.chi2a_factor=a[9];
122  node.chi2a_corr=a[10];
123  }
124 
125  // Now find the r and z boundaries that will be used during swimming
126  FindBoundaries();
127 
128  IS_VALID = true;
129 }
130 
131 //-----------------
132 // FindBoundaries
133 //-----------------
135 {
136  /// Look for boundaries where the material density changes
137  /// and record them for faster checking during swimming
138 
139  // Boundaries are recorded for r(or z) bins for which ANY adjacent
140  // z(or r) bin shows a factor of 2 or more change in the
141  // radiation length of the material.
142  //
143  // This is done by making 2 passes over every bin in the map
144  // (one for the r boundaries and one for the z) and recording
145  // the index of the bin for which an increase was detected.
146  // Doing it this way will detect a boundary multiple times
147  // (one for each bin in the direction the boundary runs).
148  // Thus, we use a set container to record the boundaries
149  // temporarily since writing the same value to it multiple
150  // times won't result in multiple entries.
151  //
152  // A final pass at the end records the r (or z) boundary values
153  // for the border between the bins in the vectors containers
154  // that are actually used by the EstimatedDistanceToBoundary
155  // method below.
156 
157  // Loop over r bins
158  set<int> iz_boundaries;
159  for(int ir=0; ir<Nr; ir++){
160 
161  // initialize boundary
162  double RadLen1 = nodes[ir][0].RadLen;
163 
164  // Loop over z bins
165  for(int iz=1; iz<Nz; iz++){
166  double RadLen2 = nodes[ir][iz].RadLen;
167  if(RadLen1<0.5*RadLen2 || RadLen2<0.5*RadLen1){
168  iz_boundaries.insert(iz);
169  }
170  }
171  }
172 
173  // Loop over z bins
174  set<int> ir_boundaries;
175  for(int iz=0; iz<Nz; iz++){
176 
177  // initialize boundary
178  double RadLen1 = nodes[0][iz].RadLen;
179 
180  // Loop over r bins
181  for(int ir=1; ir<Nr; ir++){
182  double RadLen2 = nodes[ir][iz].RadLen;
183  if(RadLen1<0.5*RadLen2 || RadLen2<0.5*RadLen1){
184  ir_boundaries.insert(ir);
185  }
186  }
187  }
188 
189  // Fill r_boundaries and z_boundaries with actual values of
190  // borders between the bins
191  for(set<int>::iterator iter=iz_boundaries.begin(); iter!=iz_boundaries.end(); iter++){
192  int iz = *iter;
193  z_boundaries.push_back(zmin+dz*(double)iz);
194  }
195  for(set<int>::iterator iter=ir_boundaries.begin(); iter!=ir_boundaries.end(); iter++){
196  int ir = *iter;
197  r_boundaries.push_back(rmin+dr*(double)ir);
198  }
199 
200 
201  // If the number of boundaries in either direction is more
202  // than half the number of bins then flag this material
203  // as having an irregular density profile.
204  irregular_density_profile = false;
205  irregular_density_profile |= (int)r_boundaries.size() > (Nr/2);
206  irregular_density_profile |= (int)z_boundaries.size() > (Nz/2);
207  if(irregular_density_profile){
208  r_boundaries.clear();
209  z_boundaries.clear();
210  }
211 }
212 
213 //-----------------
214 // FindMat
215 //-----------------
216 jerror_t DMaterialMap::FindMat(DVector3 &pos,double &rhoZ_overA, double &rhoZ_overA_logI, double &RadLen) const
217 {
218  const MaterialNode *node = FindNode(pos);
219  if(!node)return RESOURCE_UNAVAILABLE;
220 
221  rhoZ_overA = node->rhoZ_overA;
222  rhoZ_overA_logI = node->rhoZ_overA_logI;
223  RadLen = node->RadLen;
224 
225  return NOERROR;
226 }
227 
228 //-----------------
229 // FindMat
230 //-----------------
231 jerror_t DMaterialMap::FindMatALT1(DVector3 &pos,double &KrhoZ_overA,
232  double &rhoZ_overA, double &logI,
233  double &RadLen) const
234 {
235  const MaterialNode *node = FindNode(pos);
236  if(!node)return RESOURCE_UNAVAILABLE;
237 
238  KrhoZ_overA=node->KrhoZ_overA;
239  rhoZ_overA = node->rhoZ_overA;
240  logI = node->LogI;
241  RadLen = node->RadLen;
242 
243  return NOERROR;
244 }
245 
246 
247 //-----------------
248 // FindMat
249 //-----------------
250 jerror_t DMaterialMap::FindMat(DVector3 &pos, double &density, double &A, double &Z, double &RadLen) const
251 {
252  const MaterialNode *node = FindNode(pos);
253  if(!node)return RESOURCE_UNAVAILABLE;
254 
255  density = node->Density;
256  A = node->A;
257  Z = node->Z;
258  RadLen = node->RadLen;
259 
260  return NOERROR;
261 }
262 
263 // The Kalman filter needs a slightly different variant of the material
264 // properties
266  double &K_rho_Z_over_A,
267  double &rho_Z_over_A,
268  double &LogI, double &chi2c_factor,
269  double &chi2a_factor, double &chi2a_corr,
270  double &Z) const
271 {
272 
273  const MaterialNode *node = FindNode(pos);
274  if(!node)return RESOURCE_UNAVAILABLE;
275 
276  Z=node->Z;
277  rho_Z_over_A = node->rhoZ_overA;
278  LogI = node->LogI;
279  K_rho_Z_over_A=node->KrhoZ_overA;
280  chi2a_factor=node->chi2a_factor;
281  chi2a_corr=node->chi2a_corr;
282  chi2c_factor=node->chi2c_factor;
283 
284  return NOERROR;
285 }
286 
287 
288 //-----------------
289 // IsInMap
290 //-----------------
291 bool DMaterialMap::IsInMap(const DVector3 &pos) const
292 {
293  double pos_x = pos.X();
294  double pos_y = pos.Y();
295  double r = sqrt(pos_x*pos_x + pos_y*pos_y);
296  double z = pos.Z();
297  return (r>=rmin) && (r<=rmax) && (z>=zmin) && (z<=zmax);
298 }
299 
300 //-----------------
301 // EstimatedDistanceToBoundary
302 //-----------------
304 {
305  /// Give a very rough estimate of the distance a track at this position/momentum
306  /// will travel before seeing either a significant change in the raditation
307  /// length of material, or the edge of this map. The primary purpose of
308  /// this is to help determine the appropriate step size when swimming
309  /// charged tracks. As we approach a boundary of materials, we want to take
310  /// smaller steps to ensure the material is integrated properly.
311  ///
312  /// This method can be called for points either inside or outside of this map.
313  /// It can be that this map overlaps another so for points outside, we look
314  /// for the point where the track would enter this map.
315  ///
316  /// If a problem is encountered (e.g. point is outside of map and not pointing
317  /// toward it) we return a value of 1.0E6
318 
319  // The method here is to use either dr/dz or dz/dr to project a straight line
320  // in r/z space and find bins along that line that can be be checked for
321  // increases in material density. This will clearly give a worse answer than
322  // projecting a helix, but should be pretty quick and robust. Speed is critical
323  // here since this is called for every step during swimming.
324 
325  double s_to_boundary = 1.0E6;
326  if(!ENABLE_BOUNDARY_CHECK)return s_to_boundary; // low-level, catch-all opportunity to NOT do this.
327 
328  // This gets called often enough it is worth trying to bail as quickly
329  // as possible. We check if the point is outside the boundary and headed
330  // away from the map first.
331  double z = pos.Z();
332  double pz = mom.Z();
333  if(pz>0.0){
334  if(z>zmax)return s_to_boundary;
335  } else {
336  if(z<zmin)return s_to_boundary;
337  }
338 
339  double pos_x = pos.X();
340  double pos_y = pos.Y();
341  double mom_x = mom.X();
342  double mom_y = mom.Y();
343  double x_dot_p = pos_x*mom_x + pos_y*mom_y;
344  double r = sqrt(pos_x*pos_x + pos_y*pos_y);
345  if(x_dot_p>0.0){
346  if(r>rmax)return s_to_boundary;
347  } else {
348  if(r<rmin)return s_to_boundary;
349  }
350 
351  double pr = sqrt(mom_x*mom_x + mom_y*mom_y) * (x_dot_p>0 ? +1.0:-1.0); // sign of pr depends on x_dot_p
352 
353  // Unit vector pointing in momentum direction in r-z space
354  double mod = sqrt(pr*pr + pz*pz);
355  if(mod<1.0E-6)return s_to_boundary; // for when momentum is purely in phi direction
356  double p_hatR = pr/mod;
357  double p_hatZ = pz/mod;
358  //DVector2 p_hat(p_hatR, p_hatZ);
359 
360  // Get shortest distance to boundary of entire map
361  s_to_boundary = DistanceToBox(r, z, p_hatR, p_hatZ, rmin, rmax, zmin, zmax);
362 
363  // If point is outside of map, then return now. s_to_boundary has valid answer
364  if(r<=rmin || r>=rmax || z<=zmin || z>=zmax)return s_to_boundary;
365 
366  // At this point we know that the point is inside of this map. We need to
367  // search for boundaries within this map to find the closest one to the
368  // track. If the map is marked as "irregular" then we do a bin by bin
369  // search for changes in density. Otherwise, we loop over the boundary
370  // map already generated during construction.
371  if(irregular_density_profile){
372  return EstimatedDistanceToBoundarySearch(r, z, p_hatR, p_hatZ, s_to_boundary);
373  }
374 
375  // Loop over boundaries stored in z_boundaries
376  for(unsigned int i=0; i<z_boundaries.size(); i++){
377  double dist = (z_boundaries[i]-z)/p_hatZ;
378  if(!isfinite(dist))continue;
379  if(dist<0.0)continue; // boundary is behind us
380  if(dist>s_to_boundary)continue; // boundary intersection is already further than closest boundary
381 
382  // Now we must calculate whether this intersection point is
383  // still inside the map.
384  double my_r = r + dist*p_hatR;
385  if(my_r>=rmin && my_r<=rmax)s_to_boundary = dist;
386  }
387 
388  // Loop over boundaries stored in r_boundaries
389  for(unsigned int i=0; i<r_boundaries.size(); i++){
390  double dist = (r_boundaries[i]-r)/p_hatR;
391  if(!isfinite(dist))continue;
392  if(dist<0.0)continue; // boundary is behind us
393  if(dist>s_to_boundary)continue; // boundary intersection is already further than closest boundary
394 
395  // Now we must calculate whether this intersection point is
396  // still inside the map.
397  double my_z = z + dist*p_hatZ;
398  if(my_z>=zmin && my_z<=zmax)s_to_boundary = dist;
399  }
400 
401  return s_to_boundary;
402 }
403 
404 //-----------------
405 // EstimatedDistanceToBoundarySearch
406 //-----------------
407 double DMaterialMap::EstimatedDistanceToBoundarySearch(double r, double z, double p_hatR, double p_hatZ, double &s_to_boundary)
408 {
409  /// Estimate the distance to the nearest boundary (marked by a
410  /// a change in radiation length by a factor of 2 or more) by
411  /// doing an exhaustive bin-by-bin search in the direction of
412  /// the momentum vector in r-z space. This should only get called
413  /// for maps with an irregular density profile (i.e. one whose
414  /// boundaries do not appear to be in pure r or pure z directions.
415 
416  // If we got this far then the point is inside of this map. We need to check
417  // if there is a jump in the radiation length of a cell that is in the path
418  // of the trajectory so we can return that. Otherwise, we'll return the distance
419  // passed in via s_to_boundary.
420 
421  // We want to form a vector pointing in the direction of momentum in r,z space
422  // but with a magnitude such that we are taking a step across a single grid
423  // point in either r or z, whichever is smallest.
424  double scale_r=fabs(dr/p_hatR);
425  double scale_z=fabs(dz/p_hatZ);
426  double scale=1000.0;
427  if(isfinite(scale_r) && scale_r<scale)scale = scale_r;
428  if(isfinite(scale_z) && scale_z<scale)scale = scale_z;
429  DVector2 delta_rz(scale*p_hatR, scale*p_hatZ);
430 
431  // Find radiation length of our starting cell
432  int ir_start = (int)floor((r-rmin)*one_over_dr);
433  int iz_start = (int)floor((z-zmin)*one_over_dz);
434  double RadLen_start = nodes[ir_start][iz_start].RadLen;
435 
436  // Loop until we find a change of radiation length within this map or
437  // until we hit the edge of our boundaries.
438  DVector2 rzpos_start(r, z);
439  DVector2 last_rzpos = rzpos_start;
440  int last_ir = ir_start;
441  int last_iz = iz_start;
442  for(int Nsteps=0; Nsteps<MAX_BOUNDARY_SEARCH_STEPS; Nsteps++){ // limit us to looking only 10 grid points away for speed
443  // Step to next cell
444  DVector2 rzpos = last_rzpos;
445  rzpos+= delta_rz;
446 
447  // Find indexes for this cell
448  int ir = (int)floor((rzpos.X()-rmin)*one_over_dr);
449  int iz = (int)floor((rzpos.Y()-zmin)*one_over_dz);
450 
451  // Check if we hit the boundary of the map and if so, simply return
452  // the value found above for the distance to the map boundary.
453  if(ir<0 || ir>=Nr || iz<0 || iz>=Nz)return s_to_boundary;
454 
455  // Check radiation length against start point's
456  double RadLen = nodes[ir][iz].RadLen;
457  if(RadLen < 0.5*RadLen_start || RadLen>2.0*RadLen_start ){
458  double rmin_cell = (double)last_ir*dr + rmin;
459  double rmax_cell = rmin_cell + dr;
460  double zmin_cell = (double)last_iz*dz + zmin;
461  double zmax_cell = zmin_cell + dz;
462  double last_rzposX = last_rzpos.X();
463  double last_rzposY = last_rzpos.Y();
464  double s_to_cell = DistanceToBox(last_rzposX, last_rzposY, p_hatR, p_hatZ, rmin_cell, rmax_cell, zmin_cell, zmax_cell);
465  if(s_to_cell==1.0E6)s_to_cell = 0.0;
466  double my_s_to_boundary = (last_rzpos-rzpos_start).Mod() + s_to_cell;
467  if(my_s_to_boundary < s_to_boundary)s_to_boundary = my_s_to_boundary;
468 
469  return s_to_boundary;
470  }
471 
472  // Need to take another step. Update "last" variables
473  last_iz = iz;
474  last_ir = ir;
475  last_rzpos = rzpos;
476  }
477 
478  return s_to_boundary;
479 }
480 
481 //-----------------
482 // DistanceToBox
483 //-----------------
484 inline double DMaterialMap::DistanceToBox(double &posx, double &posy, double &xdir, double &ydir, double xmin, double xmax, double ymin, double ymax)
485 {
486  /// Given a point in 2-D space, a direction, and the limits of a box, find the
487  /// closest distance to box edge in the given direction.
488 
489  // Calculate intersection distances to all 4 boundaries of map
490  double dist_x1 = (xmin-posx)/xdir;
491  double dist_x2 = (xmax-posx)/xdir;
492  double dist_y1 = (ymin-posy)/ydir;
493  double dist_y2 = (ymax-posy)/ydir;
494 
495  // Make list of all positive, finite distances that are
496  // on border of box.
497  double shortestDist = 1.0E6;
498  if(isfinite(dist_x1) && dist_x1>=0.0 && dist_x1 < shortestDist){
499  //double y = (pos + dist_x1*dir).Y();
500  double y = posy + dist_x1*ydir;
501  if(y>=ymin && y<=ymax) shortestDist = dist_x1;
502  }
503  if(isfinite(dist_x2) && dist_x2>=0.0 && dist_x2 < shortestDist){
504  //double y = (pos + dist_x2*dir).Y();
505  double y = posy + dist_x2*ydir;
506  if(y>=ymin && y<=ymax) shortestDist = dist_x2;
507  }
508  if(isfinite(dist_y1) && dist_y1>=0.0 && dist_y1 < shortestDist){
509  //double x = (pos + dist_y1*dir).X();
510  double x = posx + dist_y1*xdir;
511  if(x>=xmin && x<=xmax) shortestDist = dist_y1;
512  }
513  if(isfinite(dist_y2) && dist_y2>=0.0 && dist_y2 < shortestDist){
514  //double x = (pos + dist_y2*dir).X();
515  double x = posx + dist_y2*xdir;
516  if(x>=xmin && x<=xmax) shortestDist = dist_y2;
517  }
518 
519  // Return shortest distance
520  return shortestDist;
521 }
522 
void FindBoundaries(void)
jerror_t FindMatALT1(DVector3 &pos, double &KrhoZ_overA, double &rhoZ_overA, double &logI, double &RadLen) const
double rmax
TVector2 DVector2
Definition: DVector2.h:9
TVector3 DVector3
Definition: DVector3.h:14
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
#define y
double EstimatedDistanceToBoundary(const DVector3 &pos, const DVector3 &mom)
double zmax
int Nr
Definition: bfield2root.cc:25
double rmin
int Nz
Definition: bfield2root.cc:27
double DistanceToBox(double &x, double &y, double &xdir, double &ydir, double xmin, double xmax, double ymin, double ymax)
#define _DBG_
Definition: HDEVIO.h:12
jerror_t FindMatKalman(const DVector3 &pos, double &K_rho_Z_over_A, double &rho_Z_over_A, double &LogI, double &chi2c_factor, double &chi2a_factor, double &chi2a_corr, double &Z) const
double sqrt(double)
bool IsInMap(const DVector3 &pos) const
Double_t ymin
Definition: bcal_hist_eff.C:89
Double_t ymax
Definition: bcal_hist_eff.C:91
double zmin
jerror_t FindMat(DVector3 &pos, double &rhoZ_overA, double &rhoZ_overA_logI, double &RadLen) const
double EstimatedDistanceToBoundarySearch(double r, double z, double p_hatR, double p_hatZ, double &s_to_boundary)