WireCell::Response Namespace Reference

Namespaces
	Schema

Classes
class	ColdElec
	A functional object caching gain and shape. More...
class	Generator
class	HfFilter
class	LfFilter
class	SimpleRC
class	SysResp

Functions
Schema::FieldResponse	wire_region_average (const Schema::FieldResponse &fr)
Schema::FieldResponse	average_1D (const Schema::FieldResponse &fr)
Array::array_xxf	as_array (const Schema::PlaneResponse &pr)
Array::array_xxf	as_array (const Schema::PlaneResponse &pr, int set_nrows, int set_ncols)
double	coldelec (double time, double gain=7.8, double shaping=1.0 *units::us)
	The cold electronics response function. More...
double	hf_filter (double freq, double sigma=1, double power=2, bool zero_freq_removal=true)
double	lf_filter (double freq, double tau=0.02)

Function Documentation

Array::array_xxf WireCell::Response::as_array

(

const Schema::PlaneResponse &

pr

)

Return the plane's response as a 2D array. This is a straight copy of the plane's current vectors into rows of the returned array. The first "path" will be in row 0. Each column thus contains the same sample (aka tick) across all current waveforms. Column 0 is first sample. The plane response input data is taken at face value an not attempt to resolve any implicit symmetries is made.

Definition at line 279 of file Response.cxx.

{
    int nrows = pr.paths.size();
    int ncols = pr.paths[0].current.size();
    Array::array_xxf ret(nrows, ncols); // warning, uninitialized

    for (int irow = 0; irow < nrows; ++irow) {
        auto& path = pr.paths[irow];
        for (int icol = 0; icol < ncols; ++icol) {
            ret(irow,icol) = path.current[icol]; // maybe there is a fast way to do this copy?
        }
    }
    return ret;        
}

Array::array_xxf WireCell::Response::as_array	(	const Schema::PlaneResponse &	pr,
		int	set_nrows,
		int	set_ncols
	)

Definition at line 256 of file Response.cxx.

{
    int nrows = pr.paths.size();
    int ncols = pr.paths[0].current.size();
    Array::array_xxf ret= Array::array_xxf::Zero(set_nrows, set_ncols); // warning, uninitialized

    if (set_nrows< nrows || set_ncols < ncols){
        error("Response: array dimension not correct! ");
        return ret;
    }
    
    for (int irow = 0; irow < nrows; ++irow) {
      if (irow < set_nrows){
        auto& path = pr.paths[irow];
        for (int icol = 0; icol < ncols; ++icol) {
          if (icol < set_ncols)
            ret(irow,icol) = path.current[icol]; // maybe there is a fast way to do this copy?
        }
      }
    }
    return ret;        
}

Response::Schema::FieldResponse WireCell::Response::average_1D

(

const Schema::FieldResponse &

fr

)

Definition at line 224 of file Response.cxx.

{
  using namespace WireCell::Waveform;
  using namespace WireCell::Response::Schema;

  FieldResponse fr_wire_avg = Response::wire_region_average(fr);
  
  std::vector<PlaneResponse> newplanes;
  for (auto plane : fr_wire_avg.planes) {
    std::vector<PathResponse> newpaths;

    int nsamples = Response::as_array(plane).cols();
    
    realseq_t ave_response(nsamples,0);
    
    for (auto path : plane.paths) {
      for (int k=0;k!=nsamples;k++){
        ave_response.at(k) += path.current.at(k);
      }
    }

    newpaths.push_back(PathResponse(ave_response,0.0,0.0));      
        
    newplanes.push_back(PlaneResponse(newpaths,
                                      plane.planeid,
                                      plane.location,
                                      plane.pitch));
  }
  return FieldResponse(newplanes, fr.axis, fr.origin, fr.tstart, fr.period, fr.speed);
}

double WireCell::Response::coldelec	(	double	time,
		double	gain = 7.8,
		double	shaping = 1.0*units::us
	)

The cold electronics response function.

Definition at line 350 of file Response.cxx.

{
    if (time <=0 || time >= 10 * units::microsecond) { // range of validity
        return 0.0;
    }

    const double reltime = time/shaping;

    // a scaling is needed to make the anti-Lapalace peak match the expected gain
    // fixme: this scaling is slightly dependent on shaping time.  See response.py
    gain *= 10*1.012;

    return 4.31054*exp(-2.94809*reltime)*gain
        -2.6202*exp(-2.82833*reltime)*cos(1.19361*reltime)*gain
        -2.6202*exp(-2.82833*reltime)*cos(1.19361*reltime)*cos(2.38722*reltime)*gain
        +0.464924*exp(-2.40318*reltime)*cos(2.5928*reltime)*gain
        +0.464924*exp(-2.40318*reltime)*cos(2.5928*reltime)*cos(5.18561*reltime)*gain
        +0.762456*exp(-2.82833*reltime)*sin(1.19361*reltime)*gain
        -0.762456*exp(-2.82833*reltime)*cos(2.38722*reltime)*sin(1.19361*reltime)*gain
        +0.762456*exp(-2.82833*reltime)*cos(1.19361*reltime)*sin(2.38722*reltime)*gain
        -2.620200*exp(-2.82833*reltime)*sin(1.19361*reltime)*sin(2.38722*reltime)*gain 
        -0.327684*exp(-2.40318*reltime)*sin(2.5928*reltime)*gain + 
        +0.327684*exp(-2.40318*reltime)*cos(5.18561*reltime)*sin(2.5928*reltime)*gain
        -0.327684*exp(-2.40318*reltime)*cos(2.5928*reltime)*sin(5.18561*reltime)*gain
        +0.464924*exp(-2.40318*reltime)*sin(2.5928*reltime)*sin(5.18561*reltime)*gain;
}

double WireCell::Response::hf_filter	(	double	freq,
		double	sigma = 1,
		double	power = 2,
		bool	zero_freq_removal = true
	)

Definition at line 377 of file Response.cxx.

                                                                            {
  if (flag){
    if (freq==0) return 0;
  }

  return exp(-0.5*pow(freq/sigma,power));
  
}

double WireCell::Response::lf_filter	(	double	freq,
		double	tau = 0.02
	)

Definition at line 386 of file Response.cxx.

                                                 {
  return 1-exp(-pow(freq/tau,2));
}

Response::Schema::FieldResponse WireCell::Response::wire_region_average

(

const Schema::FieldResponse &

fr

)

Return a reduced FieldResponse structure where the Path::Response::current arrays are reduced by averaging over each wire region.

Warning! this function is NOT GENERAL. It is actually specific to Garfield 1D line of paths with half the impact positions represented!

Definition at line 99 of file Response.cxx.

{
    using namespace WireCell::Waveform;
    using namespace WireCell::Response::Schema;

    std::vector<PlaneResponse> newplanes;
    for (auto plane : fr.planes) {
        std::vector<PathResponse> newpaths;

        double pitch = plane.pitch;

        std::map<int, realseq_t> avgs;
        //std::map<int, int> nums;

        
        std::map<int,realseq_t> fresp_map;
        std::map<int,std::pair<double,double>> pitch_pos_range_map;
         
        // figure out the range of each response ... 

        int nsamples=0;
        
        for (auto path : plane.paths) {
          int eff_num = path.pitchpos/(0.01 * pitch);
          if (fresp_map.find(eff_num) == fresp_map.end()){
            fresp_map[eff_num] = path.current;
          }else{
            nsamples = path.current.size();
            for (int k=0;k!=nsamples;k++){
              fresp_map[eff_num].at(k) = (fresp_map[eff_num].at(k) + path.current.at(k))/2.;
            }
          }
          if (fresp_map.find(-eff_num) == fresp_map.end()){
            fresp_map[-eff_num] = path.current;
          }else{
            int nsamples = path.current.size();
            for (int k=0;k!=nsamples;k++){
              fresp_map[-eff_num].at(k) = (fresp_map[-eff_num].at(k) + path.current.at(k))/2.;
            }
          }
        }


        std::vector<int> pitch_pos;
        for (auto it = fresp_map.begin(); it!= fresp_map.end(); it++){
          pitch_pos.push_back((*it).first);
        }
        
        double min = -1e9;
        double max = 1e9;
        for (size_t i=0;i!=pitch_pos.size();i++){
          int eff_num = pitch_pos.at(i);
          if (i==0){
            pitch_pos_range_map[eff_num] = std::make_pair(min,(pitch_pos.at(i) + pitch_pos.at(i+1))/2.*0.01*pitch);
          }else if (i==pitch_pos.size()-1){
            pitch_pos_range_map[eff_num] = std::make_pair((pitch_pos.at(i) + pitch_pos.at(i-1))/2.*0.01*pitch,max);
          }else{
            pitch_pos_range_map[eff_num] = std::make_pair((pitch_pos.at(i) + pitch_pos.at(i-1))/2.*0.01*pitch,(pitch_pos.at(i) + pitch_pos.at(i+1))/2.*0.01*pitch);
          }
        }


        // figure out how many wires ...
        std::set<int> wire_regions;
        for (size_t i=0;i!=pitch_pos.size();i++){
          if (pitch_pos.at(i)>0){
            wire_regions.insert( round((pitch_pos.at(i)*0.01*pitch-0.001*pitch)/pitch));
          }else{
            wire_regions.insert( round((pitch_pos.at(i)*0.01*pitch+0.001*pitch)/pitch));
          }
        }
        

        // do the average ... 
        for(auto it = wire_regions.begin(); it!=wire_regions.end(); it++){
          int wire_no = *it;
          if (avgs.find(wire_no) == avgs.end()) {
            avgs[wire_no] = realseq_t(nsamples);
          }
          for (auto it1 =  fresp_map.begin(); it1!= fresp_map.end(); it1++){
            int resp_num = (*it1).first;
            realseq_t& response = (*it1).second;
            double low_limit = pitch_pos_range_map[resp_num].first;
            double high_limit = pitch_pos_range_map[resp_num].second;
            if (low_limit < (wire_no - 0.5)*pitch ){
              low_limit = (wire_no - 0.5)*pitch;
            }
            if (high_limit > (wire_no+0.5)*pitch ){
              high_limit = (wire_no+0.5)*pitch;
            }

            //
            
            if (high_limit > low_limit){
              for (int k=0;k!=nsamples;k++){
                avgs[wire_no].at(k) += response.at(k) * (high_limit - low_limit) / pitch;
              }
            }
          }
        }
        
        
        // do average.
        for (auto it : avgs) {
          int region = it.first;
          realseq_t& response = it.second;

          double sum = 0;
          for (int k=0;k!=nsamples;k++){
            sum += response.at(k);
          }
          
          // pack up everything for return.
          newpaths.push_back(PathResponse(response, region*pitch, 0.0));
        }
        newplanes.push_back(PlaneResponse(newpaths,
                                          plane.planeid,
                                          plane.location,
                                          plane.pitch));
    }
    return FieldResponse(newplanes, fr.axis, fr.origin, fr.tstart, fr.period, fr.speed);
}