AttenCurve.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Store the results of the attenuation fits
3 /// \author bckhouse@caltech.edu
4 /// \date
5 /////////////////////////////////////////////////////////////////////////
6 // C/C++ includes
7 #include <cassert>
8 #include <cstdio>
9 #include <iostream>
10 #include <set>
11 
12 #include "TString.h"
13 
14 // Framework includes
16 
17 // NOvASoft includes
19 
21 
22 namespace calib
23 {
24  //......................................................................
26  {
27  AttenCurve* res = new AttenCurve;
28  res->det = det;
29  res->chan = chan;
30  res->initialized = false;
31  res->center_offset = 0;
32  return res;
33  }
34 
35  template<class T> T FOURTH(T x){return x*x*x*x;}
36  //......................................................................
37  float AttenCurve::MeanPEPerCmAt(double w) const
38  {
39  return MeanPEPerCmAt(w,"");
40  }
41 
42  //......................................................................
44  {
46 
47  //Now check just once in Calibrator instead
48  //if(!IsCalibrated()) return -1;
49 
50  // TODO: take account of center_offset here for muon catcher planes. Need
51  // to fix the todo below at the same time.
52 
53  //The attenuation formulae has a minor correction relative to tag v6 and earlier.
54  //Backwards compatibility is ensured in calibrator by passing the tag version to
55  //this function. Users can also call this function directly without specifying a
56  //tag in which case they get the latest, corrected, version of the attenuation
57  //formulae.
58 
59  // Distance to travel the long way round the fibre
60  double wlong;
61  if(version=="v1"||version=="v2"||version=="v3"||version=="v3.1"||version=="v4"||version=="v5"||version=="v5.1"||version=="v6"){
62  wlong = 1.5*cell_length+w;
63  }
64  else{
65  wlong = cell_length+w;
66  }
67 
68  const double expPart = coeff_exp*(exp(w/atten_length)+
69  exp(-wlong/atten_length))+background;
70 
71  double edgePart = 1;
72  if(w < -edge_low) edgePart -= coeff_low *FOURTH(w+edge_low);
73  if(w > edge_high) edgePart -= coeff_high*FOURTH(w-edge_high);
74 
75  double interpPart = 1;
76  if(interp_pts.size() >= 2){
77  typedef std::vector<Pt>::const_iterator it_t;
78 
79  // Find the point after w
80  it_t itNext = std::lower_bound(interp_pts.begin(), interp_pts.end(), w);
81  if(itNext == interp_pts.end()) --itNext;
82  if(itNext == interp_pts.begin()) ++itNext;
83  // And the point before
84  it_t itPrev = itNext; --itPrev;
85 
86  // Interpolate
87  const double w0 = itPrev->w; const double y0 = itPrev->factor;
88  const double w1 = itNext->w; const double y1 = itNext->factor;
89 
90  interpPart = ((w1-w)*y0+(w-w0)*y1)/(w1-w0);
91  }
92 
93  return expPart*edgePart*interpPart;
94  }
95 
96  //......................................................................
98  {
99  // This isn't actually a chisq, but more an average fractional deviation of
100  // the fit from the data. 0.2 is the eyeballed position of the tails of the
101  // distribution in FD data.
102  return chisq >= 0 && chisq < 0.2;
103 
104  // Checking the LOWESS points is another way to achieve something
105  // similar. But this implementation tends to give to much weight to
106  // problems that only occur in single locations.
107 
108  // for(unsigned int n = 0; n < interp_pts.size(); ++n){
109  // // Don't let the interpolated corrections pull more than 15% off the
110  // // original fit. 15% is arbitrary, based on a desire to calibrate a
111  // // reasonable number of channels.
112  // if(fabs(interp_pts[n].factor-1) > .15){
113  // // Don't count interpolations where just the final point is off as bad
114  // if(n > 0 && n < interp_pts.size()-1)
115  // return false;
116  // }
117  // }
118  // return true;
119  }
120 
121  //......................................................................
122  void AttenCurve::AddInterpPoint(float w, float factor)
123  {
124  Pt pt = {w, factor};
125  interp_pts.push_back(pt);
126  }
127 
128  //......................................................................
129  void AttenCurve::WriteToCSVs(FILE* fConsts, FILE* fPoints, bool mc) const
130  {
131  const int vldTime = 1; // 1970
132 
133 
135  fprintf(fConsts, "%g,%g,%g,%g,%g,%g,%g,%g,%d,%d\n",
138  chan.ToDBValidityChan(), vldTime);
139 
140  const int N = interp_pts.size();
141  for(int n = 0; n < N; ++n){
142  fprintf(fPoints, "%g,%g,%d,%d\n",
143  interp_pts[n].w, interp_pts[n].factor,
144  1000*chan.ToDBValidityChan()+n, vldTime);
145  }
146  }
147 
148  //......................................................................
149  std::ostream& operator<<(std::ostream& os, const AttenCurve& res)
150  {
151  os << "Attenuation fit results for "
153  << ", " << res.chan << ": ";
154 
155  if(!res.initialized){
156  os << " Uninitialized." << std::endl;
157  return os;
158  }
159 
160  os << res.atten_length << "cm " << res.coeff_exp << " amplitude on "
161  << res.background << " background. Rolloffs start at -"
162  << res.edge_low << "cm and +" << res.edge_high << "cm with size "
163  << res.coeff_low << " and " << res.coeff_high
164  << " fit quality " << res.chisq << std::endl;
165  os << "Interpolation points are: [ ";
166  const int N = res.interp_pts.size();
167  for(int n = 0; n < N; ++n)
168  os << "(" << res.interp_pts[n].w << ", "
169  << res.interp_pts[n].factor << ") ";
170  os << "]" << std::endl;
171  return os;
172  }
173 
174 } // end namespace calib
175 //////////////////////////////////////////////////////////////////////////
static std::string GetName(int id)
T FOURTH(T x)
Definition: AttenCurve.cxx:35
Float_t y1[n_points_granero]
Definition: compare.C:5
bool IsCalibrated() const
Definition: AttenCurve.cxx:97
static AttenCurve * Uninitialized(int det, geo::OfflineChan chan)
Return a new AttenCurve objects with fields uninitialized.
Definition: AttenCurve.cxx:25
float center_offset
Nonzero in short muon-catcher cells. Positive is closer to readout.
Definition: AttenCurve.h:51
int ToDBValidityChan() const
Definition: OfflineChan.cxx:19
CDPStorage service.
Encapsulate the geometry of one entire detector (near, far, ndos)
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
const XML_Char * version
Definition: expat.h:187
geo::OfflineChan chan
Definition: AttenCurve.h:53
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A (plane, cell) pair.
Definition: OfflineChan.h:17
assert(nhit_max >=nhit_nbins)
float MeanPEPerCmAt(double w) const
Mean response of this channel at this distance from detector centre.
Definition: AttenCurve.cxx:37
double T
Definition: Xdiff_gwt.C:5
Float_t w
Definition: plot.C:20
void WriteToCSVs(FILE *fConsts, FILE *fPoints, bool mc) const
Definition: AttenCurve.cxx:129
std::ostream & operator<<(std::ostream &os, const AttenCurve &res)
Definition: AttenCurve.cxx:149
std::vector< Pt > interp_pts
Definition: AttenCurve.h:43
void AddInterpPoint(float w, float factor)
Definition: AttenCurve.cxx:122