Lots of new stuff

doc/wup.tex

@@ -952,4 +952,241 @@ M_{00}^{-1}
 +(2M_{01}M_{12}M_{02} -M_{00}M_{12}^2 -M_{11}M_{02}^2)/M_{22}}
 \end{eqnarray}
 
+\section{2D interpolation}
+Using eight nearest neighbours around a central bin, then the
+coefficients of a Taylor expansion to second order (which has six
+coefficients) can be found using a fit. Specifically
+\begin{eqnarray}
+  E(x,y)
+  &=& E(0,0)
+  + \left.\frac{\partial E}{\partial x}\right|_{00} x
+  + \left.\frac{\partial E}{\partial y}\right|_{00} y
+  + \left.\frac{\partial^2 E}{\partial x^2}\right|_{00} \frac{x^2}{2}
+  + \left.\frac{\partial^2 E}{\partial x\partial y}\right|_{00} xy
+  + \left.\frac{\partial^2 E}{\partial y^2}\right|_{00} \frac{y^2}{2}\\
+  &=& E_0 + E_x x + E_y y + E_{xx} x^2 + E_{xy} xy + E_{yy} y^2
+\end{eqnarray}
+The values in the eight bins at $x$ and $y$ within $\pm 1$ of the central
+bin are
+\begin{eqnarray}
+  e(-1, 1) = e_1 &=& E_0 - E_x + E_y + E_{xx} - E_{xy} + E_{yy}\\
+  e( 0, 1) = e_2 &=& E_0 + E_y + E_{yy}\\
+  e( 1, 1) = e_3 &=& E_0 + E_x + E_y + E_{xx} + E_{xy} + E_{yy}\\
+  e(-1, 0) = e_4 &=& E_0 - E_x + E_{xx}\\
+  e( 1, 0) = e_5 &=& E_0 + E_x + E_{xx}\\
+  e(-1,-1) = e_6 &=& E_0 - E_x - E_y + E_{xx} + E_{xy} + E_{yy}\\
+  e( 0,-1) = e_7 &=& E_0 - E_y + E_{yy}\\
+  e( 1,-1) = e_8 &=& E_0 + E_x - E_y + E_{xx} - E_{xy} + E_{yy}
+\end{eqnarray}
+so
+\begin{eqnarray}
+\begin{pmatrix}
+e_1 \\ e_2 \\ e_3 \\ e_4 \\ e_5 \\ e_6 \\ e_7 \\ e_8
+\end{pmatrix}
+= \begin{pmatrix}
+ \phantom{-}1 & -1 &\phantom{-}1 &\phantom{-}1 & -1 &\phantom{-}1\\
+ \phantom{-}1 &\phantom{-}0 &\phantom{-}1 &\phantom{-}0 &\phantom{-}0 &\phantom{-}1\\
+ \phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1\\
+ \phantom{-}1 & -1 &\phantom{-}0 &\phantom{-}1 &\phantom{-}0 &\phantom{-}0\\
+ \phantom{-}1 & \phantom{-}1 &\phantom{-}0 &\phantom{-}1 &\phantom{-}0 &\phantom{-}0\\
+ \phantom{-}1 & -1 & -1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1\\
+ \phantom{-}1 & \phantom{-}0 & -1 &\phantom{-}0 &\phantom{-}0 & \phantom{-}1\\
+ \phantom{-}1 & \phantom{-}1 & -1 &\phantom{-}1 & -1 &\phantom{-}1 
+\end{pmatrix}
+\begin{pmatrix}
+E_0 \\ E_x \\ E_y \\ E_{xx} \\ E_{xy} \\ E_{yy}
+\end{pmatrix}
+\end{eqnarray}
+or $\underline{e}=\underline{\underline{a}}\underline{E}$. The problem is then
+to find $\underline{E}$ given $\underline{e}$. This is an overconstrained
+problem so a fit is required. Forming
+\begin{eqnarray}
+\chi^2 = \sum_{k=1}^8 \left(e_k - \sum_{j=1}^6 a_{kj}E_j\right)^2
+\end{eqnarray}
+then
+\begin{eqnarray}
+  \frac{\partial \chi^2}{\partial E_i} = \sum_{k=1}^8
+  - 2\left(e_k - \sum_{j=1}^6 a_{kj} E_j \right) a_{ki}
+\end{eqnarray}
+and so for the minimum $\chi^2$
+\begin{eqnarray}
+  \sum_{j=1}^6  \sum_{k=1}^8 a_{ki} a_{kj} E_j
+  = \sum_{k=1}^8 e_k a_{ki}
+  \qquad\mbox{or}\qquad
+  \underline{\underline{a}}^T \underline{\underline{a}} \underline{E}
+  = \underline{\underline{a}}^T \underline{e}
+\end{eqnarray}
+With $\underline{\underline{M}}=\underline{\underline{a}}^T \underline{\underline{a}}$, then
+\begin{eqnarray}
+  \underline{E}
+  = \underline{\underline{M}}^{-1} \underline{\underline{a}}^T \underline{e}
+\end{eqnarray}
+gives the required solution. Using $\underline{\underline{a}}$ above, then
+\begin{eqnarray}
+  \underline{\underline{M}}
+  &=& \begin{pmatrix}
+ \phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1 & \phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1\\
+ -1 &\phantom{-}0 &\phantom{-}1 & -1 &\phantom{-}1 & -1 &\phantom{-}0 &\phantom{-}1\\
+ \phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}0 &\phantom{-}0 & -1 & -1 & -1 \\
+ \phantom{-}1 & \phantom{-}0 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1 & \phantom{-}1 & \phantom{-}0 &\phantom{-}1\\
+ -1 &\phantom{-}0 &\phantom{-}1 & \phantom{-}0 &\phantom{-}0 & \phantom{-}1 &\phantom{-}0 & -1\\
+ \phantom{-}1 &\phantom{-}1 &\phantom{-}1 & \phantom{-}0 &\phantom{-}0 &\phantom{-}1 & \phantom{-}1 &\phantom{-}1
+  \end{pmatrix}
+  \begin{pmatrix}
+ \phantom{-}1 & -1 &\phantom{-}1 &\phantom{-}1 & -1 &\phantom{-}1\\
+ \phantom{-}1 &\phantom{-}0 &\phantom{-}1 &\phantom{-}0 &\phantom{-}0 &\phantom{-}1\\
+ \phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1\\
+ \phantom{-}1 & -1 &\phantom{-}0 &\phantom{-}1 &\phantom{-}0 &\phantom{-}0\\
+ \phantom{-}1 & \phantom{-}1 &\phantom{-}0 &\phantom{-}1 &\phantom{-}0 &\phantom{-}0\\
+ \phantom{-}1 & -1 & -1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1\\
+ \phantom{-}1 & \phantom{-}0 & -1 &\phantom{-}0 &\phantom{-}0 & \phantom{-}1\\
+ \phantom{-}1 & \phantom{-}1 & -1 &\phantom{-}1 & -1 &\phantom{-}1
+  \end{pmatrix}\\
+  &=&   \begin{pmatrix}
+ \phantom{-}8 &\phantom{-}0 &\phantom{-}0 &\phantom{-}6 &\phantom{-}0 &\phantom{-}6\\
+ \phantom{-}0 &\phantom{-}6 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0\\
+ \phantom{-}0 &\phantom{-}0 &\phantom{-}6 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0\\
+ \phantom{-}6 &\phantom{-}0 &\phantom{-}0 &\phantom{-}6 &\phantom{-}0 &\phantom{-}4\\
+ \phantom{-}0 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0 &\phantom{-}4 &\phantom{-}0\\
+ \phantom{-}6 &\phantom{-}0 &\phantom{-}0 &\phantom{-}4 &\phantom{-}0 &\phantom{-}6
+  \end{pmatrix}
+\end{eqnarray}
+for which
+\begin{eqnarray}
+  \underline{\underline{M}}^{-1}
+  &=&   \begin{pmatrix}
+ \phantom{-}5/4 &\phantom{-}0 &\phantom{-}0 & -3/4 &\phantom{-}0 & -3/4\\
+ \phantom{-}0 &\phantom{-}1/6 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0\\
+ \phantom{-}0 &\phantom{-}0 &\phantom{-}1/6 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0\\
+ -3/4 &\phantom{-}0 &\phantom{-}0 &\phantom{-}3/4 &\phantom{-}0 &\phantom{-}1/4\\
+ \phantom{-}0 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0 &\phantom{-}1/4 &\phantom{-}0\\
+ -3/4 &\phantom{-}0 &\phantom{-}0 &\phantom{-}1/4 &\phantom{-}0 &\phantom{-}3/4
+  \end{pmatrix}
+\end{eqnarray}
+Hence
+\begin{eqnarray}
+  \underline{\underline{M}}^{-1}\underline{\underline{a}}^T
+  &=&   \begin{pmatrix}
+ \phantom{-}5/4 &\phantom{-}0 &\phantom{-}0 & -3/4 &\phantom{-}0 & -3/4\\
+ \phantom{-}0 &\phantom{-}1/6 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0\\
+ \phantom{-}0 &\phantom{-}0 &\phantom{-}1/6 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0\\
+ -3/4 &\phantom{-}0 &\phantom{-}0 &\phantom{-}3/4 &\phantom{-}0 &\phantom{-}1/4\\
+ \phantom{-}0 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0 &\phantom{-}1/4 &\phantom{-}0\\
+ -3/4 &\phantom{-}0 &\phantom{-}0 &\phantom{-}1/4 &\phantom{-}0 &\phantom{-}3/4
+  \end{pmatrix}
+\begin{pmatrix}
+ \phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1 & \phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1\\
+ -1 &\phantom{-}0 &\phantom{-}1 & -1 &\phantom{-}1 & -1 &\phantom{-}0 &\phantom{-}1\\
+ \phantom{-}1 &\phantom{-}1 &\phantom{-}1 &\phantom{-}0 &\phantom{-}0 & -1 & -1 & -1 \\
+ \phantom{-}1 & \phantom{-}0 &\phantom{-}1 &\phantom{-}1 &\phantom{-}1 & \phantom{-}1 & \phantom{-}0 &\phantom{-}1\\
+ -1 &\phantom{-}0 &\phantom{-}1 & \phantom{-}0 &\phantom{-}0 & \phantom{-}1 &\phantom{-}0 & -1\\
+ \phantom{-}1 &\phantom{-}1 &\phantom{-}1 & \phantom{-}0 &\phantom{-}0 &\phantom{-}1 & \phantom{-}1 &\phantom{-}1
+\end{pmatrix}\nonumber\\
+&=& \begin{pmatrix}
+ -1/4 &\phantom{-}1/2 & -1/4 & \phantom{-}1/2 & \phantom{-}1/2 & -1/4 & \phantom{-}1/2 & -1/4\\
+ -1/6 &\phantom{-}0 & \phantom{-}1/6 & -1/6 & \phantom{-}1/6 & -1/6 & \phantom{-}0 & \phantom{-}1/6\\
+ \phantom{-}1/6 &\phantom{-}1/6 & \phantom{-}1/6 &\phantom{-}0 & \phantom{-}0 & -1/6 & -1/6 & -1/6\\
+ \phantom{-}1/4 & -1/2 & \phantom{-}1/4  &\phantom{-}0 & \phantom{-}0 & \phantom{-}1/4 & -1/2 & \phantom{-}1/4\\
+ -1/4 & \phantom{-}0 & \phantom{-}1/4  &\phantom{-}0 & \phantom{-}0 & \phantom{-}1/4 & \phantom{-}0 & -1/4\\
+ \phantom{-}1/4 & \phantom{-}0 & \phantom{-}1/4  & -1/2 & -1/2 & \phantom{-}1/4 & \phantom{-}0 & \phantom{-}1/4
+  \end{pmatrix}
+\end{eqnarray}
+Hence, the interpolation to the central bin is simply done by the first line of the above, i.e.
+\begin{eqnarray}
+E_0 = \frac{e_2+e_4+e_5+e_7}{2} - \frac{e_1+e_3+e_6+e_8}{4}
+\end{eqnarray}
+where the first group are the horizonal and vertical nearest neighbours and
+the second group are the diagonal nearest neighbours.
+
+\section{Trigger rates}
+Consider a target fraction of BXs which trigger of $1/r$, so $r$ is the
+average number of BX between triggers. The fat BX period is $f$
+and the time multiplex period is $m$. For current assumptions,
+$1/r = 750/40000 = 0.01875$ so $r = 53.33$, while $f=107$ and $m=18$.
+
+With a fat event prescale $s$, then for finite $s$ the target fraction
+must be maintained. Let $p=1/a$ be the probability of a non-fat BX triggering.
+Within one fat BX period, there is one fat BX and $f-1$ non-fat BXs. The fat BX
+can trigger due to the usual L1Accept or due to prescaling. The total
+probability of it triggering is
+\begin{eqnarray}
+1 - (1- p)\left(1- \frac{1}{s}\right)
+= p + \frac{1}{s} - \frac{p}{s}
+\end{eqnarray}
+Hence the fraction of BXs which trigger is
+\begin{eqnarray}
+p\frac{(f-1)}{f} + \left(p + \frac{1}{s} - \frac{p}{s}\right)\frac{1}{f}
+= p + \frac{1-p}{sf}
+\end{eqnarray}
+which must be $1/r$. Hence
+\begin{eqnarray}
+\frac{f}{r} = p\left[f - \frac{1}{s}\right]+ \frac{1}{s}
+\end{eqnarray}
+and so 
+\begin{eqnarray}
+p = \frac{1}{a} = \frac{f/r - 1/s}{f - 1/s}
+\end{eqnarray}
+Clearly for $s = \infty$ then $p = 1/a = 1/r$ as expected. Note, this requires
+$f > r$ for values of the prescale down to $s=1$; i.e. the fat BX period must
+not be less than the average number of BXs between triggers for a prescale of
+one. Note also that the above is independent of any time multiplexing.
+
+Considering time multiplexing, the assuming the fat BX period is a prime, then
+the repetition period is $mf$ BXs long. Within this time, there are $m$ fat
+BXs and $f$ time multiplex BXs. Of these, one BX is both fat and time
+multiplexed, so there are $m-1$ fat BXs which are not time multiplexed and
+$f-1$ time multiplexed BXs which are not fat. The remainder, which is
+\begin{eqnarray}
+mf -(m-1)-(f-1)-1 = mf-m-f+1 = m(f-1)-(f-1) = (m-1)(f-1)
+\end{eqnarray}
+are normal BXs. The number of the $mf$ BXs which trigger is
+\begin{eqnarray}
+(p + 1/s - p/s)\frac{1}{mf},\quad
+  (p + 1/s - p/s)\frac{m-1}{mf},\quad
+  p\frac{f-1}{mf},\quad
+  p\frac{(m-1)(f-1)}{mf}
+\end{eqnarray}
+On a TPG time multiplexing board which handles all L1Accepts but only saves
+non-trivial data for time multiplexed BXs, then the second and fourth
+categories are treated as normal L1Accepts, with a total fraction of
+\begin{eqnarray}
+F_0 &=& (p + 1/s - p/s)\frac{m-1}{mf}+ p\frac{(m-1)(f-1)}{mf}
+= p \frac{(m-1)+(m-1)(f-1)}{mf} + \frac{(1-p)}{s}\frac{(m-1)}{mf}\\
+&=& p \frac{(m-1)f}{mf} + \frac{(1-p)}{s}\frac{(m-1)}{mf}
+= p \frac{m-1}{m} + \frac{(1-p)}{s}\frac{(m-1)}{mf}
+\end{eqnarray}
+The time multiplexed events have a fraction
+\begin{eqnarray}
+F_1 =  p\frac{f-1}{mf}
+\end{eqnarray}
+and the fat time multiplexed events have a fraction
+\begin{eqnarray}
+F_2 = (p + 1/s - p/s)\frac{1}{mf}
+\end{eqnarray}
+To check, the total fraction is
+\begin{eqnarray}
+  F_0+F_1+F_2 &=& p \frac{m-1}{m} + \frac{(1-p)}{s}\frac{(m-1)}{mf}
+  +p\frac{f-1}{mf}
+  +(p + 1/s - p/s)\frac{1}{mf}\\
+&=& \frac{1}{mf}\left[pf(m-1)+ \frac{(1-p)(m-1)}{s} + p(f-1) 
+    + p + \frac{1}{s} - \frac{p}{s}\right]\\
+  &=& \frac{1}{mf}\left[pfm 
+ + \frac{1}{s}\left(1+(1-p)(m-1) - p\right)\right]\\
+&=& p + \frac{m(1-p)}{smf}
+= p + \frac{(1-p)}{sf}
+\end{eqnarray}
+as shown above. For the assumption numbers gives previously and $s=1$, then
+\begin{eqnarray}
+p = 0.009493,\quad
+  F_0 = 0.01771\quad
+  F_1 = 0.0005225,\quad
+  F_2 = 0.0005192
+\end{eqnarray}
+which in terms of rate in kHz are
+\begin{eqnarray}
+  R_0 = 708.333,\quad
+  R_1 = 20.898,\quad
+  R_2 = 20.768
+\end{eqnarray}
+
 \end{document}

interface/Coefficients.hh

@@ -199,6 +199,12 @@ public:
     return (fErg/fNum)-sum*fCoe;
   }
 
+  double energy(const double *p) {
+    TVectorD v(kNumberOfLayers);
+    for(unsigned i(0);i<kNumberOfLayers;i++) v(i)=p[i];
+    return energy(v);
+  }
+
   double energy(const TVectorD &p) const {
     assert(p.GetNrows()==kNumberOfLayers);
 

interface/Particle.hh

@@ -94,6 +94,12 @@ public:
 
     fFourP=TLorentzVector( p4.x(), p4.y(), p4.z(), p4.e() );
 
+    if(fPdgId==22 || abs(fPdgId)==12 || abs(fPdgId)==14 || abs(fPdgId)==16) fMass=0.0;
+    else {
+      fMass=p4.m();
+      //if(fMass<0.0005)
+      //std::cout << "PDG = " << fPdgId << ", mass = " << fMass << std::endl;
+    }
   }
 
   Particle(int pdg, double px, double py, double pz) : pGenParticle(0) {
@@ -113,6 +119,10 @@ public:
     return fPdgId;
   }
 
+  int particleCharge() const {
+    return charge(fPdgId);
+  }
+
   Point position() const {
     return fPosition;
   }
@@ -121,6 +131,10 @@ public:
     return fMomentum;
   }
 
+  double energy() const {
+    return sqrt(fMomentum.rSquared()+fMass*fMass);
+  }
+  
   TLorentzVector FourP() const {
     return fFourP;
   } 
@@ -137,6 +151,50 @@ public:
     return fTrgC3d;
   }
 
+  Particle projectBarrel(double b) const {
+    BACKTRACE_AND_ASSERT(b>0.0);
+
+    Particle p;
+
+    int q(charge(fPdgId));
+
+    if(q==0) {
+      double r02(fPosition.rhoSquared());
+      BACKTRACE_AND_ASSERT(b*b>r02);
+      double r0DotP(fPosition.x()*fMomentum.x()+fPosition.y()*fMomentum.y());
+      double l((-r0DotP+sqrt(r0DotP+fMomentum.rhoSquared()*(b*b-r02)))/fMomentum.rhoSquared());
+      
+      p.pGenParticle=pGenParticle;
+      p.fPdgId=fPdgId;
+      p.fPosition=fPosition+(fMomentum*l);
+      p.fMomentum=fMomentum;
+      p.fFourP=fFourP;
+
+    } else {
+      double phiLimit(acos(-1.0)*q*fBField*b/fMomentum.rho());
+      Point2D rxy;
+      for(unsigned i(0);i<100;i++) {
+	double dPhi(0.01*(i+1)*phiLimit);
+        rxy.x(fPosition.x()+(fMomentum.x()*sin(dPhi)+fMomentum.y()*(1.0-cos(dPhi)))/(q*fBField));
+        rxy.y(fPosition.y()+(fMomentum.y()*sin(dPhi)-fMomentum.x()*(1.0-cos(dPhi)))/(q*fBField));
+	if(rxy.rho()>b) {
+	  p.pGenParticle=pGenParticle;
+	  p.fPdgId=fPdgId;
+	  p.fPosition.setXY(rxy.x(),rxy.y());
+	  p.fPosition.z(fPosition.z()+fMomentum.z()*dPhi/(q*fBField));
+	  double px(fMomentum.x()*cos(dPhi)+fMomentum.y()*sin(dPhi));
+	  double py(fMomentum.y()*cos(dPhi)-fMomentum.x()*sin(dPhi));
+	  p.fMomentum=Point(px,py,fMomentum.z());
+	  p.fFourP.SetX(p.fMomentum.x());
+	  p.fFourP.SetY(p.fMomentum.y());
+	  p.fFourP.SetZ(p.fMomentum.z());
+	  return p;
+	}
+      }	
+    }
+    return p;
+  }
+
   Particle project(double z) const {
     Particle p;
     p.pGenParticle=pGenParticle;
@@ -162,8 +220,8 @@ public:
       */
 
       p.fMomentum=fMomentum;
-
       p.fFourP=fFourP;
+
     } else {
         double dPhi(dz*q*fBField/fMomentum.z());
 
@@ -202,7 +260,8 @@ private:
   int fPdgId;
   Point fPosition;
   Point fMomentum;
-
+  double fMass;
+  
   TLorentzVector fFourP; 
 
   TrgC3d *fTrgC3d;

interface/Point.hh

@@ -4,7 +4,7 @@
 #include <iostream>
 #include <cstdint>
 
-#include "Point2D.hh"
+#include "Scale2D.hh"
 
 class Point : public Point2D {
 public:
@@ -27,8 +27,12 @@ public:
     return fZ;
   }
 
+  double rSquared() const {
+    return fX*fX+fY*fY+fZ*fZ;
+  }
+
   double r() const {
-    return sqrt(fX*fX+fY*fY+fZ*fZ);
+    return sqrt(rSquared());
   }
 
   double cosTheta() const {
@@ -47,8 +51,9 @@ public:
     return asinh(z()/rho());
   }
 
-  Point2D scale() const {
-    return Point2D(fX/fZ,fY/fZ);
+  Point2D scale(bool absZ=false) const {
+    if(absZ) return Scale2D(fX/fabs(fZ),fY/fabs(fZ));
+    return Scale2D(fX/fZ,fY/fZ);
   }
 
   void flipEndcap() {
@@ -56,7 +61,11 @@ public:
     fZ=-fZ;
   }
 
-  double deltaR(const Point &p) {
+  void z(double d) {
+    fZ=d;
+  }
+  
+  double deltaR(const Point &p) const {
     double dEta(eta()-p.eta());
 
     double dPhi(phi()-p.phi());
@@ -67,7 +76,7 @@ public:
     return sqrt(dEta*dEta+dPhi*dPhi);
   }
 
-  Point operator+(const Point2D &xy) {
+  Point operator+(const Point2D &xy) const {
     return Point(fX+xy.x(),fY+xy.y(),fZ);
   } 
 
@@ -76,7 +85,7 @@ public:
     fY+=xy.y();
   } 
 
-  Point operator+(const Point &p) {
+  Point operator+(const Point &p) const {
     return Point(fX+p.x(),fY+p.y(),fZ+p.fZ);
   } 
 
@@ -86,7 +95,7 @@ public:
     fZ+=p.z();
   } 
 
-  Point operator-(const Point2D &xy) {
+  Point operator-(const Point2D &xy) const {
     return Point(fX-xy.x(),fY-xy.y(),fZ);
   } 
 
@@ -95,7 +104,7 @@ public:
     fY-=xy.y();
   } 
 
-  Point operator-(const Point &p) {
+  Point operator-(const Point &p) const {
     return Point(fX-p.x(),fY-p.y(),fZ-p.fZ);
   } 
 
@@ -105,7 +114,7 @@ public:
     fZ-=p.z();
   } 
 
-  Point operator*(const double &d) {
+  Point operator*(const double &d) const {
     return Point(fX*d,fY*d,fZ*d);
   } 
 
@@ -115,7 +124,7 @@ public:
     fZ*=d;
   } 
 
-  Point operator/(const double &d) {
+  Point operator/(const double &d) const {
     return Point(fX/d,fY/d,fZ/d);
   } 
 

interface/Point2D.hh

@@ -5,6 +5,8 @@
 #include <cstdint>
 #include <cmath>
 
+#include "Backtrace.hh"
+
 class Point2D {
 public:
   Point2D()
@@ -17,6 +19,11 @@ public:
 
   ~Point2D() {
   }
+
+  void zero() {
+    fX=0.0;
+    fY=0.0;
+  }
   
   void x(double d) {
     fX=d;
@@ -36,6 +43,12 @@ public:
     fY=r*sin(f);
   }
 
+  void setRS(double r, double s) {
+    BACKTRACE_AND_ASSERT(r>0.0);
+    fX=r*cos(s/r);
+    fY=r*sin(s/r);
+  }
+
   double x() const {
     return fX;
   }
@@ -44,8 +57,12 @@ public:
     return fY;
   }
 
+  double rhoSquared() const {
+    return fX*fX+fY*fY;
+  }
+
   double rho() const {
-    return sqrt(fX*fX+fY*fY);
+    return sqrt(rhoSquared());
   }
 
   double phi(bool positiveRange=false) const {
@@ -54,6 +71,10 @@ public:
     return f;
   }
 
+  double s() const {
+    return rhoPhi();
+  }
+
   double rhoPhi() const {
     return rho()*phi();
   }

interface/Scale2D.hh

+#ifndef Scale2D_HH
+#define Scale2D_HH
+
+#include <iostream>
+#include <cstdint>
+
+#include "Point2D.hh"
+
+class Scale2D : public Point2D {
+public:
+  Scale2D()
+    : Point2D() {
+  }
+
+  Scale2D(double x, double y)
+    : Point2D(x,y) {
+  }
+
+  Scale2D(const Point2D &xy)
+    : Point2D(xy) {
+  }
+
+  ~Scale2D() {
+  }
+  
+  double r() const {
+    return sqrt(rhoSquared()+1.0);
+  }
+
+  double cosTheta() const {
+    return 1.0/r();
+  }
+
+  double sinTheta() const {
+    return rho()/r();
+  }
+
+  double tanTheta() const {
+    return rho();
+  }
+
+  double eta() const {
+    return asinh(1.0/rho());
+  }
+
+  double deltaR(const Scale2D &p) {
+    double dEta(eta()-p.eta());
+
+    double dPhi(phi()-p.phi());
+    double pi(acos(-1.0));
+    if     (dPhi<-pi) dPhi+=2.0*pi;
+    else if(dPhi> pi) dPhi-=2.0*pi;
+
+    return sqrt(dEta*dEta+dPhi*dPhi);
+  }
+
+  Scale2D operator+(const Scale2D &xy) const {
+    return Scale2D(fX+xy.x(),fY+xy.y());
+  } 
+
+  void operator+=(const Scale2D &xy) {
+    fX+=xy.x();
+    fY+=xy.y();
+  } 
+
+  Scale2D operator-(const Scale2D &xy) const {
+    return Scale2D(fX-xy.x(),fY-xy.y());
+  } 
+
+  void operator-=(const Scale2D &xy) {
+    fX-=xy.x();
+    fY-=xy.y();
+  } 
+
+  Scale2D operator*(const double &d) const {
+    return Scale2D(fX*d,fY*d);
+  } 
+
+  void operator*=(const double &d) {
+    fX*=d;
+    fY*=d;
+  }