25 MAX_BOUNDARY_SEARCH_STEPS = 30;
26 ENABLE_BOUNDARY_CHECK =
true;
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.");
31 this->namepath = namepath;
36 this->jcalib = jcalib;
42 vector< vector<float> > Mmap;
43 jcalib->Get(namepath, Mmap);
57 map<float, int> rvals;
58 map<float, int> zvals;
62 for(
unsigned int i=0; i<Mmap.size(); i++){
63 vector<float> &a = Mmap[i];
78 dr = (rmax-
rmin)/(
double)(
Nr-1);
79 dz = (zmax-
zmin)/(
double)(
Nz-1);
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;
98 for(
int ir=0; ir<
Nr; ir++){
103 for(
unsigned int i=0; i<Mmap.size(); i++){
104 vector<float> &a = Mmap[i];
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];
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];
158 set<int> iz_boundaries;
159 for(
int ir=0; ir<
Nr; ir++){
162 double RadLen1 = nodes[ir][0].RadLen;
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);
174 set<int> ir_boundaries;
175 for(
int iz=0; iz<
Nz; iz++){
178 double RadLen1 = nodes[0][iz].RadLen;
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);
191 for(set<int>::iterator iter=iz_boundaries.begin(); iter!=iz_boundaries.end(); iter++){
193 z_boundaries.push_back(zmin+dz*(
double)iz);
195 for(set<int>::iterator iter=ir_boundaries.begin(); iter!=ir_boundaries.end(); iter++){
197 r_boundaries.push_back(rmin+dr*(
double)ir);
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();
219 if(!node)
return RESOURCE_UNAVAILABLE;
232 double &rhoZ_overA,
double &logI,
233 double &RadLen)
const
236 if(!node)
return RESOURCE_UNAVAILABLE;
253 if(!node)
return RESOURCE_UNAVAILABLE;
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,
274 if(!node)
return RESOURCE_UNAVAILABLE;
293 double pos_x = pos.X();
294 double pos_y = pos.Y();
295 double r =
sqrt(pos_x*pos_x + pos_y*pos_y);
297 return (r>=rmin) && (r<=rmax) && (z>=
zmin) && (z<=zmax);
325 double s_to_boundary = 1.0E6;
326 if(!ENABLE_BOUNDARY_CHECK)
return s_to_boundary;
334 if(z>zmax)
return s_to_boundary;
336 if(z<zmin)
return s_to_boundary;
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);
346 if(r>rmax)
return s_to_boundary;
348 if(r<rmin)
return s_to_boundary;
351 double pr =
sqrt(mom_x*mom_x + mom_y*mom_y) * (x_dot_p>0 ? +1.0:-1.0);
354 double mod =
sqrt(pr*pr + pz*pz);
355 if(mod<1.0E-6)
return s_to_boundary;
356 double p_hatR = pr/mod;
357 double p_hatZ = pz/mod;
361 s_to_boundary = DistanceToBox(r, z, p_hatR, p_hatZ, rmin, rmax, zmin, zmax);
364 if(r<=rmin || r>=rmax || z<=zmin || z>=zmax)
return s_to_boundary;
371 if(irregular_density_profile){
372 return EstimatedDistanceToBoundarySearch(r, z, p_hatR, p_hatZ, s_to_boundary);
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;
380 if(dist>s_to_boundary)
continue;
384 double my_r = r + dist*p_hatR;
385 if(my_r>=rmin && my_r<=rmax)s_to_boundary = dist;
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;
393 if(dist>s_to_boundary)
continue;
397 double my_z = z + dist*p_hatZ;
398 if(my_z>=zmin && my_z<=zmax)s_to_boundary = dist;
401 return s_to_boundary;
424 double scale_r=fabs(dr/p_hatR);
425 double scale_z=fabs(dz/p_hatZ);
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);
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;
440 int last_ir = ir_start;
441 int last_iz = iz_start;
442 for(
int Nsteps=0; Nsteps<MAX_BOUNDARY_SEARCH_STEPS; Nsteps++){
448 int ir = (int)floor((rzpos.X()-
rmin)*one_over_dr);
449 int iz = (int)floor((rzpos.Y()-
zmin)*one_over_dz);
453 if(ir<0 || ir>=Nr || iz<0 || iz>=
Nz)
return s_to_boundary;
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;
469 return s_to_boundary;
478 return s_to_boundary;
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;
497 double shortestDist = 1.0E6;
498 if(isfinite(dist_x1) && dist_x1>=0.0 && dist_x1 < shortestDist){
500 double y = posy + dist_x1*ydir;
501 if(y>=ymin && y<=ymax) shortestDist = dist_x1;
503 if(isfinite(dist_x2) && dist_x2>=0.0 && dist_x2 < shortestDist){
505 double y = posy + dist_x2*ydir;
506 if(y>=ymin && y<=ymax) shortestDist = dist_x2;
508 if(isfinite(dist_y1) && dist_y1>=0.0 && dist_y1 < shortestDist){
510 double x = posx + dist_y1*xdir;
511 if(x>=xmin && x<=xmax) shortestDist = dist_y1;
513 if(isfinite(dist_y2) && dist_y2>=0.0 && dist_y2 < shortestDist){
515 double x = posx + dist_y2*xdir;
516 if(x>=xmin && x<=xmax) shortestDist = dist_y2;
void FindBoundaries(void)
jerror_t FindMatALT1(DVector3 &pos, double &KrhoZ_overA, double &rhoZ_overA, double &logI, double &RadLen) const
double EstimatedDistanceToBoundary(const DVector3 &pos, const DVector3 &mom)
double DistanceToBox(double &x, double &y, double &xdir, double &ydir, double xmin, double xmax, double ymin, double ymax)
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
bool IsInMap(const DVector3 &pos) const
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)