#include #include #include using namespace std; #include #include #include "DMaterialMap.h" #include "DMagneticFieldMap.h" //----------------- // DMaterialMap (Constructor) //----------------- DMaterialMap::DMaterialMap(string namepath, JCalibration *jcalib) { /// Read the specified material map in from the calibration database. /// This will read in the map and figure out the number of grid /// points in each direction (r, and z) and the range in each. MAX_BOUNDARY_SEARCH_STEPS = 30; ENABLE_BOUNDARY_CHECK = true; 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(...)"); gPARMS->SetDefaultParameter("GEOM:ENABLE_BOUNDARY_CHECK", ENABLE_BOUNDARY_CHECK, "Enable boundary checking (superceeds any setting in DReferenceTrajectory). This is for debugging only."); this->namepath = namepath; // Read in map from calibration database. This should really be // built into a run-dependent geometry framework, but for now // we do it this way. this->jcalib = jcalib; if(!jcalib)return; string blank(' ',80); //cout< > Mmap; jcalib->Get(namepath, Mmap); //cout<<(int)Mmap.size()<<" entries ("; if(Mmap.size()<1){ //cout<<")"< to make a // histogram of the entries by using the key to hold the extent // so that the number of entries will be equal to the number of // different values. map rvals; map zvals; double rmin, zmin, rmax, zmax; rmin = zmin = 1.0E6; rmax = zmax = -1.0E6; for(unsigned int i=0; i &a = Mmap[i]; float &r = a[0]; float &z = a[1]; rvals[r] = 1; zvals[z] = 1; if(rrmax)rmax=r; if(z>zmax)zmax=z; } Nr = rvals.size(); Nz = zvals.size(); r0 = rmin; z0 = zmin; dr = (rmax-rmin)/(double)(Nr-1); dz = (zmax-zmin)/(double)(Nz-1); one_over_dr=1./dr; one_over_dz=1./dz; //cout<<"Nr="<rmin = rmin-dr/2.0; this->rmax = rmax+dr/2.0; this->zmin = zmin-dz/2.0; this->zmax = zmax+dz/2.0; // Set sizes of nested vectors to hold node data nodes.resize(Nr); for(int ir=0; ir &a = Mmap[i]; float &r = a[0]; float &z = a[1]; int ir = (int)floor((r-this->rmin)/dr); int iz = (int)floor((z-this->zmin)/dz); if(ir<0 || ir>=Nr){_DBG_<<"ir out of range: ir="<s_to_boundary)continue; // boundary intersection is already further than closest boundary // Now we must calculate whether this intersection point is // still inside the map. double my_r = r + dist*p_hatR; if(my_r>=rmin && my_r<=rmax)s_to_boundary = dist; } // Loop over boundaries stored in r_boundaries for(unsigned int i=0; is_to_boundary)continue; // boundary intersection is already further than closest boundary // Now we must calculate whether this intersection point is // still inside the map. double my_z = z + dist*p_hatZ; if(my_z>=zmin && my_z<=zmax)s_to_boundary = dist; } return s_to_boundary; } //----------------- // EstimatedDistanceToBoundarySearch //----------------- double DMaterialMap::EstimatedDistanceToBoundarySearch(double r, double z, double p_hatR, double p_hatZ, double &s_to_boundary) { /// Estimate the distance to the nearest boundary (marked by a /// a change in radiation length by a factor of 2 or more) by /// doing an exhaustive bin-by-bin search in the direction of /// the momentum vector in r-z space. This should only get called /// for maps with an irregular density profile (i.e. one whose /// boundaries do not appear to be in pure r or pure z directions. // If we got this far then the point is inside of this map. We need to check // if there is a jump in the radiation length of a cell that is in the path // of the trajectory so we can return that. Otherwise, we'll return the distance // passed in via s_to_boundary. // We want to form a vector pointing in the direction of momentum in r,z space // but with a magnitude such that we are taking a step across a single grid // point in either r or z, whichever is smallest. double scale_r=fabs(dr/p_hatR); double scale_z=fabs(dz/p_hatZ); double scale=1000.0; if(isfinite(scale_r) && scale_r=Nr || iz<0 || iz>=Nz)return s_to_boundary; // Check radiation length against start point's double RadLen = nodes[ir][iz].RadLen; if(RadLen < 0.5*RadLen_start || RadLen>2.0*RadLen_start ){ double rmin_cell = (double)last_ir*dr + rmin; double rmax_cell = rmin_cell + dr; double zmin_cell = (double)last_iz*dz + zmin; double zmax_cell = zmin_cell + dz; double last_rzposX = last_rzpos.X(); double last_rzposY = last_rzpos.Y(); double s_to_cell = DistanceToBox(last_rzposX, last_rzposY, p_hatR, p_hatZ, rmin_cell, rmax_cell, zmin_cell, zmax_cell); if(s_to_cell==1.0E6)s_to_cell = 0.0; double my_s_to_boundary = (last_rzpos-rzpos_start).Mod() + s_to_cell; if(my_s_to_boundary < s_to_boundary)s_to_boundary = my_s_to_boundary; return s_to_boundary; } // Need to take another step. Update "last" variables last_iz = iz; last_ir = ir; last_rzpos = rzpos; } return s_to_boundary; } //----------------- // DistanceToBox //----------------- inline double DMaterialMap::DistanceToBox(double &posx, double &posy, double &xdir, double &ydir, double xmin, double xmax, double ymin, double ymax) { /// Given a point in 2-D space, a direction, and the limits of a box, find the /// closest distance to box edge in the given direction. // Calculate intersection distances to all 4 boundaries of map double dist_x1 = (xmin-posx)/xdir; double dist_x2 = (xmax-posx)/xdir; double dist_y1 = (ymin-posy)/ydir; double dist_y2 = (ymax-posy)/ydir; // Make list of all positive, finite distances that are // on border of box. double shortestDist = 1.0E6; if(isfinite(dist_x1) && dist_x1>=0.0 && dist_x1 < shortestDist){ //double y = (pos + dist_x1*dir).Y(); double y = posy + dist_x1*ydir; if(y>=ymin && y<=ymax) shortestDist = dist_x1; } if(isfinite(dist_x2) && dist_x2>=0.0 && dist_x2 < shortestDist){ //double y = (pos + dist_x2*dir).Y(); double y = posy + dist_x2*ydir; if(y>=ymin && y<=ymax) shortestDist = dist_x2; } if(isfinite(dist_y1) && dist_y1>=0.0 && dist_y1 < shortestDist){ //double x = (pos + dist_y1*dir).X(); double x = posx + dist_y1*xdir; if(x>=xmin && x<=xmax) shortestDist = dist_y1; } if(isfinite(dist_y2) && dist_y2>=0.0 && dist_y2 < shortestDist){ //double x = (pos + dist_y2*dir).X(); double x = posx + dist_y2*xdir; if(x>=xmin && x<=xmax) shortestDist = dist_y2; } // Return shortest distance return shortestDist; }