Various studies

doc/wup.tex

@@ -537,6 +537,38 @@ Representation & Exponent & Mantissa & Value \cr\hline
 \end{tabular}
 \end{center}
 
+To get an unbiased decoded value, it is necessary to decode to a value a factor
+of two larger than the original value. To do this, all values stored with the
+exponent greater than 1 have 0b01111\dots appended to the mantissa where the
+number of 1 bits is set to be one less than the original LSB. Expanding the
+example table above, this gives
+
+\bigskip
+\begin{center}
+\begin{tabular}{c|c|c|c|c|c|c|c}
+\hline
+Representation & Exponent & Mantissa & Value & $2\times$Rep & Appended & $2\times$Value & Value\cr\hline
+ 0 & 0b00 & 0b00 & 0b00000 = \phantom{2}0 & \phantom{2}0 &  & \phantom{2}0 & 0 \cr
+ 1 & 0b00 & 0b01 & 0b00001 = \phantom{2}1 & \phantom{2}2 & & \phantom{2}2 & 1 \cr
+ 2 & 0b00 & 0b10 & 0b00010 = \phantom{2}2 & \phantom{2}4 & & \phantom{2}4 & 2\cr
+ 3 & 0b00 & 0b11 & 0b00011 = \phantom{2}3 & \phantom{2}6 & & \phantom{2}6 & 3\cr
+ 4 & 0b01 & 0b00 & 0b00100 = \phantom{2}4 & \phantom{2}8 & & \phantom{2}8 & 4\cr
+ 5 & 0b01 & 0b01 & 0b00101 = \phantom{2}5 &           10 & &           10 & 5\cr
+ 6 & 0b01 & 0b10 & 0b00110 = \phantom{2}6 &           12 & &           12 & 6\cr
+ 7 & 0b01 & 0b11 & 0b00111 = \phantom{2}7 &           14 & &           14 & 7 \cr
+ 8 & 0b10 & 0b00 & 0b01000 = \phantom{2}8 &           16 & 0b0001 & 0b010001 = 17 & 8.5 \cr
+ 9 & 0b10 & 0b01 & 0b01010 = 10 &           20 & 0b0101 & 0b010101 = 21 & 10.5 \cr
+10 & 0b10 & 0b10 & 0b01100 = 12 &           24 & 0b1001 & 0b011001 = 25 & 12.5\cr
+11 & 0b10 & 0b11 & 0b01110 = 14 &           28 & 0b1101 & 0b011101 = 29 & 14.5 \cr
+12 & 0b11 & 0b00 & 0b10000 = 16 &           32 & 0b00011 & 0b0100011 = 35 & 17.5\cr
+13 & 0b11 & 0b01 & 0b10100 = 20 &           40 & 0b01011 & 0b0101011 = 43 & 21.5\cr
+14 & 0b11 & 0b10 & 0b11000 = 24 &           48 & 0b10011 & 0b0110011 = 51 & 25.5\cr
+15 & 0b11 & 0b11 & 0b11100 = 28 &           56 & 0b11011 & 0b0111011 = 59 & 29.5\cr
+\hline
+\end{tabular}
+\end{center}
+
+
 \section{Template fit of energy in depth}
 Assume a template shape for a photon of $P_l$ per photon energy GeV 
 for layer $l$.
@@ -731,6 +763,10 @@ a factor of $\sin\theta$, which is given by
 \begin{equation} 
   \sin\theta = \sin\left(\tan^{-1}\rho_s\right)
   = \sin\left[\tan^{-1}\left(\frac{1}{\sinh\eta}\right)\right]
+\qquad\mbox{Note also}\qquad
+\sin\theta = \frac{1}{\cosh\eta},\quad
+\tan\theta = \frac{1}{\sinh\eta},\quad
+\cos\theta = \tanh\eta
 \end{equation}
 %
 \begin{figure}[ht!]
@@ -1189,4 +1225,130 @@ which in terms of rate in kHz are
   R_2 = 20.768
 \end{eqnarray}
 
+\section{Photon model}
+For photons at shower maximum (around layer 10, counting from 0), the
+cumulative 0.683
+and 0.900 radii are around 1.2\,cm and 2.8\,cm, respectively.
+
+\section{Gaussian}
+Assume a Gaussian radial distribution for the shower shape in a given
+layer with a width parameter $a$, so having a density
+\begin{eqnarray}
+\rho = \frac{1}{2\pi a^2} e^{-r^2/2a^2}\,rdr\,d\phi
+\end{eqnarray}
+Check the normalisation: integrating to some radius $R$ gives
+\begin{eqnarray}
+  \int_0^R \int_0^{2\pi} \rho\,rdr\,d\phi
+  = \int_0^R \frac{r}{a^2} e^{-r^2/2a^2}\,dr
+  = \left[-e^{-r^2/2a^2}\right]_0^R = 1 - e^{-R^2/2a^2}
+\end{eqnarray}
+This goes to 1 as $R\rightarrow\infty$ and so is normalised.
+For a cumulative amount $C<1$, then the radius needed is the solution of
+\begin{eqnarray}
+  e^{-R^2/2a^2} = 1-C,
+  \qquad\mbox{so}\qquad
+  \frac{R}{a} = \sqrt{-2\ln(1-C)}
+\end{eqnarray}
+For the standard values,
+$C=0.683$ requires $R/a = 1.52$ while $C=0.900$ requires $R/a = 2.15$.
+These imply
+parameter values of $a = 1.2/1.52 = 0.79$\,cm and $a = 2.8/2.15 = 1.30$\,cm,
+respectively. Hence, they are not really compatible, but an average value of
+1.05\,cm is probably adequate.
+%
+\begin{figure}[ht!]
+\begin{center}
+\includegraphics[width=9cm]{GaussianCumulative.pdf}
+\label{fig:ExponentialCumulative}
+\end{center}
+\end{figure}
+
+\section{Exponential}
+Assume an exponential radial distribution for the shower shape in a given
+layer with an exponential parameter $a$, so having a density
+\begin{eqnarray}
+\rho = \frac{1}{2\pi a^2} e^{-r/a}\,rdr\,d\phi
+\end{eqnarray}
+Check the normalisation: integrating to some radius $R$ gives
+\begin{eqnarray}
+  \int_0^R \int_0^{2\pi} \rho\,rdr\,d\phi
+  = \int_0^R \frac{r}{a^2} e^{-r/a}\,dr
+\end{eqnarray}
+Consider
+\begin{eqnarray}
+  \frac{d}{dr}\left(\frac{r}{a}e^{-r/a}\right)
+  = \frac{1}{a}e^{-r/a} - \frac{r}{a^2}e^{-r/a}
+\end{eqnarray}
+Hence
+\begin{eqnarray}
+\int_0^R \frac{r}{a^2}e^{-r/a}\,dr
+= \int_0^R \frac{1}{a}e^{-r/a}\,dr - \left[\frac{r}{a}e^{-r/a}\right]_0^R
+= \left[-e^{-r/a}\right]_0^R - \frac{R}{a}e^{-R/a}
+= 1 - \left(1 + \frac{R}{a}\right) e^{-R/a}
+\end{eqnarray}
+This goes to 1 as $R\rightarrow\infty$ and so is normalised.
+For a cumulative amount $C<1$, then the radius needed is the solution of
+\begin{eqnarray}
+\left(1 + \frac{R}{a}\right) e^{-R/a} = 1-C
+\end{eqnarray}
+which can be solved numerically. For the standard values,
+$C=0.683$ requires $R/a = 2.36$ while $C=0.900$ requires $R/a = 3.89$.
+These imply
+parameter values of $a = 1.2/2.36 = 0.508$\,cm and $a = 2.8/3.89 = 0.720$\,cm,
+respectively. Hence, they are not really compatible, but an average value of
+0.61\,cm is probably adequate.
+%
+\begin{figure}[ht!]
+\begin{center}
+\includegraphics[width=9cm]{ExponentialCumulative.pdf}
+\label{fig:ExponentialCumulative}
+\end{center}
+\end{figure}
+
+\section{Exponential square-root}
+Assume a radial distribution of an exponential with a square root
+for the shower shape in a given
+layer with a exponential parameter $\sqrt{a}$, so having a density
+\begin{eqnarray}
+\rho = \frac{1}{24\pi a^2} e^{-\sqrt{r/a}}\,rdr\,d\phi
+\end{eqnarray}
+Check the normalisation: integrating to some radius $R$ gives
+\begin{eqnarray}
+  \int_0^R \int_0^{2\pi} \rho\,rdr\,d\phi
+  = \int_0^R \frac{r}{12a^2} e^{-\sqrt{r/a}}\,dr
+\end{eqnarray}
+Let $x=r/a$ and consider
+\begin{eqnarray}
+  \frac{d}{dx}\left[(2x^{3/2}+6x+12x^{1/2}+12)e^{-\sqrt{x}}\right]
+  &=& (3x^{1/2}+6+6x^{-1/2})e^{-\sqrt{x}}
+  -(x^{3/2}+3x+6x^{1/2}+6)x^{-1/2} e^{-\sqrt{x}}\\
+  &=& (3x^{1/2}+6+6x^{-1/2} - x-3x^{1/2}-6-6x^{-1/2})e^{-\sqrt{x}}\\
+  &=& -x e^{-\sqrt{x}}
+\end{eqnarray}
+Hence
+\begin{eqnarray}
+\int_0^R \frac{r}{12a^2} e^{-\sqrt{r/a}}\,dr
+= \left[-\left(1+x^{1/2}+\frac{x}{2}+\frac{x^{3/2}}{6}\right)e^{-\sqrt{x}}\right]_0^{R/a}
+= 1 - \left[1+(R/a)^{1/2}+\frac{R/a}{2}+\frac{(R/a)^{3/2}}{6} \right] e^{-\sqrt{R/a}}
+\end{eqnarray}
+This goes to 1 as $R\rightarrow\infty$ and so is normalised.
+For a cumulative amount $C<1$, then the radius needed is the solution of
+\begin{eqnarray}
+\left[1+(R/a)^{1/2}+\frac{R/a}{2}+\frac{(R/a)^{3/2}}{6} \right] e^{-\sqrt{R/a}}
+ = 1-C
+\end{eqnarray}
+which can be solved numerically. For the standard values,
+$C=0.683$ requires $R/a = 21.7$ while $C=0.900$ requires $R/a = 44.6$.
+These imply
+parameter values of $a = 1.2/21.7 = 0.0553$\,cm and $a = 2.8/44.6 = 0.0628$\,cm,
+respectively. Hence, they are quite compatible, and an average value of
+0.059\,cm is probably adequate.
+%
+\begin{figure}[ht!]
+\begin{center}
+\includegraphics[width=9cm]{ExponentialSqrtCumulative.pdf}
+\label{fig:ExponentialSqrtCumulative}
+\end{center}
+\end{figure}
+
 \end{document}

interface/PackerBase.hh

@@ -8,24 +8,27 @@
 
 #include "utilities.cc"
 
-class PackerBase {
+template <class Uint>
+class PackerBaseUint {
 public:
-  PackerBase() : fDatum(0) {
+  PackerBaseUint() : fDatum(0) {
   }
  
-  ~PackerBase() {
+  ~PackerBaseUint() {
   }
   
-  void datum(uint8_t d) {
+  void datum(Uint d) {
     fDatum=d;
   }
 
-  uint8_t datum() const {
+  Uint datum() const {
     return fDatum;
   }
 
   virtual std::string name() const {
-    return "PackerBase";
+    std::ostringstream sout;
+    sout << "PackerBaseUint" << sizeof(Uint);
+    return sout.str();
   }
   
   virtual void pack(uint32_t d) {
@@ -36,14 +39,16 @@ public:
     assert(false);
   }
   void print() const {
-    std::cout << "PackerBase::print()  Datum = "
+    std::cout << "PackerBaseUint::print()  Datum = "
 	      << printDec(fDatum) << " = "
 	      << printHex(fDatum) << " = "
 	      << printBin(fDatum) << std::endl;
   }
 
 protected:
-  uint8_t fDatum;
+  Uint fDatum;
 };
 
+typedef PackerBaseUint<uint8_t> PackerBase;
+
 #endif

interface/PackerFloat22.hh

+#ifndef PackerFloat22_HH
+#define PackerFloat22_HH
+
+#include <iostream>
+#include <cstdint>
+#include <cmath>
+#include <cassert>
+
+#include "PackerBase.hh"
+
+class PackerFloat22 : public PackerBaseUint<uint16_t> {
+public:
+  PackerFloat22(unsigned e, unsigned m) : fMan(m), fExp(e), fBit(fMan+(1<<fExp)-1), fDataMax((uint64_t(1)<<fBit)-1) {
+    assert(m<=8*sizeof(uint16_t));
+    
+    std::cout << name() << "::ctor() Bits in mantissa = " << fMan
+	      << ", exponent = " << fExp << ", total range = "
+	      << fBit << ", data maximum = 0x" << std::hex << fDataMax
+	      << std::dec << " = " << fDataMax << std::endl;
+  }
+ 
+  ~PackerFloat22() {
+  }
+  
+  virtual std::string name() const {
+    std::ostringstream sout;
+    sout << "PackerFloat22E" << fExp << "M" << fMan;
+    return sout.str();
+  }
+  
+  void pack(uint32_t d) {
+    uint32_t dSave(d);
+
+    assert(d<=fDataMax);
+
+    if(d<(uint32_t(1)<<fMan)) {
+      fDatum=d;
+      return;
+    }
+
+    for(unsigned i(1);i<(uint32_t(1)<<fExp);i++) {
+      if(d<(uint32_t(1)<<(fMan+1))) {
+	fDatum=d-(1<<fMan)+(i<<fMan);
+	return;
+      }
+      d>>=1;
+    }
+
+    std::cout << name() << " ERROR packing value = " 
+	      << printDec(dSave) << " = "
+	      << printHex(dSave) << " = "
+	      << printBin(dSave) << std::endl;
+    assert(false);
+    return;
+  }
+
+  uint32_t unpack() const {
+    unsigned e(fDatum>>fMan);
+    assert(e<(uint32_t(1)<<fExp));
+    unsigned m(fDatum-(e<<fMan));
+    assert(m<(uint32_t(1)<<fMan));
+
+    if(e==0) return m;
+    m+=(1<<fMan);
+    if(e==1) return m;
+    m=(m<<1)+1;
+    return (m<<(e-2));
+  }
+
+  uint32_t unpackTimes2() const {
+    unsigned e(fDatum>>fMan);
+    assert(e<(uint32_t(1)<<fExp));
+    unsigned m(fDatum-(e<<fMan));
+    assert(m<(uint32_t(1)<<fMan));
+
+    if(e==0) return 2*m;
+    m+=(1<<fMan);
+    if(e==1) return 2*m;
+
+    m=(m<<16)+0x7fff;
+    return (m>>(16-e));    
+  }
+
+  uint32_t dataMax() const {
+    return fDataMax;
+  }
+  
+  void print() const {
+    std::cout << "PackerFloat22::print()";PackerBaseUint<uint16_t>::print();
+    uint32_t up(unpack());
+    std::cout << " Unpack = " 
+	      << printDec(up) << " = "
+	      << printHex(up) << " = "
+	      << printBin(up) << std::endl;
+  }
+
+protected:
+  const unsigned fMan;
+  const unsigned fExp;
+  const unsigned fBit;
+  const uint32_t fDataMax;
+};
+
+#endif

interface/TE1F.hh

@@ -9,7 +9,8 @@
 
 class TE1F {
 public:
-  TE1F(const std::string &l, const std::string &t, unsigned n, double lo, double hi) {
+  TE1F(const std::string &l, const std::string &t, unsigned n, double lo, double hi, bool u=true) {
+    up=u;
     wSum=0.0;
     r=new TH1F((l+"Raw").c_str(),t.c_str(),n,lo,hi);
     i=new TH1F((l+"Int").c_str(),t.c_str(),n,lo,hi);
@@ -23,10 +24,17 @@ public:
     r->Fill(x,w);
 
     int n(i->FindBin(x));
-    for(int b(1);b<n;b++) {
-      i->Fill(i->GetBinCenter(b),w);
-    }
-    if(x>i->GetBinCenter(n)) i->Fill(i->GetBinCenter(n),w);
+    if(up) {
+      for(int b(1);b<n;b++) {
+	i->Fill(i->GetBinCenter(b),w);
+      }
+      if(x>i->GetBinCenter(n)) i->Fill(i->GetBinCenter(n),w);
+    } else {
+      for(int b(n+1);b<=i->GetNbinsX();b++) {
+	i->Fill(i->GetBinCenter(b),w);
+      }
+      if(x<i->GetBinCenter(n)) i->Fill(i->GetBinCenter(n),w);
+    }      
 
     wSum+=w;
     for(int b(1);b<=i->GetNbinsX();b++) {
@@ -35,6 +43,7 @@ public:
   }
 
 protected:
+  bool up;
   double wSum;
   TH1F *r;
   TH1F *i;

src/DaqDaqBuffer.cpp

+/*
+
+  TO DO:
+  - Add latency to lpGBT (and other?) to get throttle timing correct
+  - Delay ECON RX->main buffer to be just in time so as to minimise buffer occ
+  - Add option for L1Accept every 44 BX: permanent continuous
+ */
 #include "Hack.hh"
 
 #include <iostream>