ScatteringSurfaces.cxx
Go to the documentation of this file.
1 ///
2 /// \file ScatteringSurfaces.cxx
3 /// \brief Construct a set of scattering surfaces for Lutz.
4 /// \author messier@indiana.edu
5 ///
7 #include <iostream>
8 #include <cmath>
11 using namespace bpfit;
12 
14  double dstol,
15  double dxtol,
16  double dttol) :
17  fPDG(abs(PDG)),
18  fdStol(dstol),
19  fdXtol(dxtol),
20  fdTtol(dttol)
21 { }
22 
23 //......................................................................
24 
26  const double* p0,
27  const double* p1,
28  const std::vector<double>& s)
29 {
30  fS. clear();
31  fSigAlpha.clear();
32 
33  Path path(geo, p0, p1, s[0]);
34  this->MakeSteps(path, s);
35 
36  if (fS.size()==0) this->ForceOnePlane(path, s);
37 }
38 
39 //......................................................................
40 ///
41 /// \brief Make one scattering surface between first and last
42 /// measurements
43 ///
44 /// \param path - The track data
45 /// \param s - Measurement plane locations
46 ///
48  const std::vector<double>& s)
49 {
50  unsigned int i1=0;
51  unsigned int i2=s.size()-1;
52  double ds, smid, dx, psi, dt;
53  this->GetStepData(path, s, i1, i2, &ds, &smid, &dx, &psi, &dt);
54  fS. push_back(smid);
55  fSigAlpha.push_back(psi);
56 }
57 
58 //......................................................................
59 ///
60 /// \brief Check the step data to see if its within tolerance
61 ///
62 /// \param ds - step size in cm
63 /// \param dx - step size in radiation lengths
64 /// \param dt - RMS transverse scattering distance in cm
65 /// \return true=good step, false=step out of tolerance.
66 ///
67 bool ScatteringSurfaces::GoodStep(double ds, double dx, double dt)
68 {
69  if (ds>fdStol) return false;
70  if (dx>fdXtol) return false;
71  if (dt>fdTtol) return false;
72  return true;
73 }
74 
75 //......................................................................
76 ///
77 /// \brief Round off near track end can cause the momentum to go to
78 /// exactly zero or even a little negative. Protect against that.
79 ///
80 /// beta - Input beta value. Possibly modified on return
81 /// p - Input momentum value. Possibly modified on return
82 ///
83 void ScatteringSurfaces::CheckLimits(double* beta, double* p)
84 {
85  static const double betatol = 1.0E-6;
86  static const double ptol = 1.0E-5;
87  if (*p<=ptol) {
88  if (fabs(*p)<ptol) *p = ptol;
89  else
90  throw BPException(__FILE__, __LINE__, kBAD_P, *p);
91  }
92  if (*beta<betatol) {
93  if (fabs(*beta)<betatol) *beta = betatol;
94  else
95  throw BPException(__FILE__, __LINE__, kBAD_P, *beta);
96  }
97 }
98 
99 //......................................................................
100 ///
101 /// \brief Get relevant data about one interval between measurements
102 ///
103 /// \param path - The track parameters
104 /// \param s - Locations of measurement planes (cm)
105 /// \param i1 - step starts at s[i1]
106 /// \param i2 - step ends at s[i2]
107 /// \param ds - on return, path length between i2 and i2 in cm
108 /// \param dx - on return, ds in units of radiation lengths
109 /// \param psi - RMS in-plane scattering angle (radians)
110 /// \param dt - RMS in-place scattering distance (cm)
111 ///
113  const std::vector<double>& s,
114  unsigned int i1,
115  unsigned int i2,
116  double* ds,
117  double* smid,
118  double* dx,
119  double* psi,
120  double* dt)
121 {
122  static const double psi0 = 13.6E-3;
123  double s1 = s[i1];
124  double s2 = s[i2];
125 
126  *smid = 0.5*(s1+s2);
127 
128  double beta=0;
129  double p =0;
130  if (fPDG==13) path.MuonParams (*smid, 0, &p, &beta, 0);
131  else if (fPDG==211) path.PionParams (*smid, 0, &p, &beta, 0);
132  else if (fPDG==2212) path.ProtonParams(*smid, 0, &p, &beta, 0);
133  else {
134  std::cout << "\n\n\nERROR: BreakPointFitter is unable to fit under the assumption PDG = "
135  << fPDG
136  << ". Aborting...\n\n\n";
137  abort();
138  }
139  this->CheckLimits(&beta, &p);
140 
141  *ds = s2-s1;
142  *dx = path.X0(s2)-path.X0(s1);
143 
144  if (*dx>1.0E-6) *psi = (psi0/beta/p)*sqrt(*dx)*(1.0+0.038*log(*dx));
145  else{
146  *psi = 1.0E-6;
147  }
148  *dt = (*psi)*(*ds);
149 
150  if(*psi < 0.0){
151  std::cout<<"\n beta: "<<beta
152  <<"\n p: "<<p
153  <<"\n dx: "<<*dx
154  <<std::endl;
155  }
156 
157 }
158 
159 //......................................................................
160 ///
161 /// \brief Make the largest step that is within tolerance
162 ///
163 /// \param path - Track parameters
164 /// \param s - Measurement plane locations
165 /// \param i1 - step starts at s[i1]
166 /// \return i2 - end point of step if s[i2]
167 ///
169  const std::vector<double>& s,
170  unsigned int i1)
171 {
172  unsigned int i2; // Index of end of step
173  double ds; // Step size in cm
174  double smid; // Midpoint of step
175  double psi; // RMS in-plane scattering angle
176  double dx; // Step size in radiation lengths
177  double dt; // Predicted RMS transverse scatter for step
178  //
179  // Increase i2 until the step data goes out of tolerance
180  //
181  for (i2=i1+1;i2+1<s.size();++i2) {
182  this->GetStepData(path, s, i1, i2, &ds, &smid, &dx, &psi, &dt);
183  if (this->GoodStep(ds, dx, dt)==false) {
184  //
185  // Last step took us out of tolerance, back up if we can
186  //
187  if (i2-i1>1) --i2;
188  break;
189  }
190  }
191  return i2;
192 }
193 
194 //......................................................................
195 ///
196 /// \brief Find good locations for the scattering surfaces
197 ///
198 /// \param path - Track parameters
199 /// \param s - Measurement plane locations (cm)
200 ///
202  const std::vector<double>& s)
203 {
204  unsigned int i1;
205  unsigned int i2;
206  double ds; // Step size in cm
207  double smid; // Path length at middle of step
208  double dx; // Step size in radiation lengths
209  double psi; // RMS in-plane scattering angle
210  double dt; // RMS transvers scattering at end of step in cm
211  for (i1=0;(i1+1)<s.size();i1=i2) {
212  i2 = this->MakeStep(path, s, i1);
213  this->GetStepData(path, s, i1, i2, &ds, &smid, &dx, &psi, &dt);
214  if(ds==0) continue;
215 
216  fS. push_back(smid);
217  fSigAlpha.push_back(psi);
218 
219  }
220 }
221 
222 //......................................................................
223 
224 unsigned int ScatteringSurfaces::N() const { return fS.size(); }
225 
226 //......................................................................
227 
228 double ScatteringSurfaces::S(unsigned int i) const { return fS[i]; }
229 
230 //......................................................................
231 
232 double ScatteringSurfaces::SigAlpha(unsigned int i) const
233 {
234  return fSigAlpha[i];
235 }
236 
237 ////////////////////////////////////////////////////////////////////////
double X0(double s) const
Return the number of radiation lengths crossed by the path.
Definition: Path.cxx:112
double fdStol
Largest allowed step (cm)
unsigned int MakeStep(const Path &path, const std::vector< double > &s, unsigned int i1)
Make the largest step that is within tolerance.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const unsigned int kBAD_P
negative value for p
Definition: BPException.h:25
void MakeSteps(const Path &path, const std::vector< double > &s)
Find good locations for the scattering surfaces.
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
vector< vector< double > > clear
bool GoodStep(double ds, double dx, double dt)
Check the step data to see if its within tolerance.
void abs(TH1 *hist)
void ProtonParams(double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
Definition: Path.cxx:462
Double_t beta
unsigned int N() const
The number of scattering planes constructed.
const XML_Char * s
Definition: expat.h:262
void MakeSurfaces(const geo::GeometryBase &geo, const double *p0, const double *p1, const std::vector< double > &s)
Build scattering surfaces.
Definition: Cand.cxx:23
void GetStepData(const Path &path, const std::vector< double > &s, unsigned int i1, unsigned int i2, double *ds, double *smid, double *dx, double *psi, double *dt)
Get relevant data about one interval between measurements.
double dx[NP][NC]
void CheckLimits(double *beta, double *p)
Round off near track end can cause the momentum to go to exactly zero or even a little negative...
void PionParams(double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
Definition: Path.cxx:449
ScatteringSurfaces(int PDG, double dstol, double dxtol, double dttol)
Construct tool for building scattering surfaces.
Track parameters (s, X0, KE, beta, ...) along a particular path in the detector.
Definition: Path.h:18
Float_t E
Definition: plot.C:20
base_types push_back(int_type())
std::vector< double > fSigAlpha
RMS scattering uncertainties.
A very simple service to remember what detector we&#39;re working in.
Construct a set of scattering surfaces for Lutz.
void MuonParams(double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
Return the estimated muon parameters at path length s.
Definition: Path.cxx:436
OStream cout
Definition: OStream.cxx:6
const std::string path
Definition: plot_BEN.C:43
unsigned int fPDG
Assumed track PDG code.
double fdXtol
Largest allowed step (rad. len.)
std::vector< double > fS
Locations of scattering surfaces.
The geometry of one entire detector (near, far, ipnd)
Definition: GeometryBase.h:49
double fdTtol
Largest allowed RMS scatter (cm)
Helper for AttenCurve.
Definition: Path.h:10
void ForceOnePlane(const Path &path, const std::vector< double > &s)
Make one scattering surface between first and last measurements.
double SigAlpha(unsigned int i) const
The predicted RMS scattering angle at plane i.
double S(unsigned int i) const
The location of the ith scattering plane.