Skip to content

GitLab

Commit 4a65be59 authored by dauncey's avatar dauncey
Browse files
Options

Many changes from the last four months

parent e74839d6
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 2006 additions and 128 deletions
doc/Float.xlsx 0 → 100644
View file @ 4a65be59
File added
doc/wup.pdf
View file @ 4a65be59
No preview for this file type
doc/wup.tex
View file @ 4a65be59
......@@ -174,9 +174,12 @@ where the positive sign is required for $q>100$\,fC.
\section{Link data representation}
The selected trigger cell
data on the link need to be represented in a small number of bits $n$,
data are calculated to a large number of bits, typically 16-18.
On the links, they need to be represented in a small number of bits $n$,
typically $\sim 8$.
This could be linear or logarithmic.
This could be linear or logarithmic or floating.
\subsection{Linear representation}
For linear, then in general it can be linear betwen $x_\mathrm{min}$ and
$x_\mathrm{max}$ and 0 or $2^n-1$ outside this range. This can be represented
within the range by
......@@ -187,8 +190,10 @@ which can be inverted to give
\begin{equation}
x = x_\mathrm{min} + \frac{y(x_\mathrm{max}-x_\mathrm{min})}{2^n-1}
\end{equation}
\subsection{Logarithmic representation}
For logarithmic, the general case would be $y=a\log(x)+b$ but with
$c=b/a$ and $x_\mathrm{min}=e^-c$, then
$c=b/a$ and $x_\mathrm{min}=e^{-c}$, then
\begin{equation}
y = a\log(x)+b = a\log(x)+ac = a(\log(x)+c) = a(\log(x)-\log(x_\mathrm{min}))
= a \log(x/x_\mathrm{min})
......@@ -202,6 +207,142 @@ This can be inverted to give
x = x_\mathrm{min}\left(\frac{x_\mathrm{max}}{x_\mathrm{min}}\right)^{y/(2^n-1)}
\end{equation}
\subsection{Float representation}
Here, the $2^n$ values are split into an exponent of $E$ bits and a mantissa of
$M$ bits. The naive approach is simply to take the actual value as the
mantissa shifted up by $E$ bits. For example, for $E=2$ and $M=2$, then the
16 possible values would give the table below.
\bigskip
\begin{center}
\begin{tabular}{c|c|c|c}
\hline
Representation & Exponent & Mantissa & Value \cr\hline
0 & 0b00 & 0b00 & 0b00000 = \phantom{2}0\cr
1 & 0b00 & 0b01 & 0b00001 = \phantom{2}1\cr
2 & 0b00 & 0b10 & 0b00010 = \phantom{2}2\cr
3 & 0b00 & 0b11 & 0b00011 = \phantom{2}3\cr
4 & 0b01 & 0b00 & 0b00000 = \phantom{2}0\cr
5 & 0b01 & 0b01 & 0b00100 = \phantom{2}2\cr
6 & 0b01 & 0b10 & 0b01000 = \phantom{2}4\cr
7 & 0b01 & 0b11 & 0b01100 = \phantom{2}6\cr
8 & 0b10 & 0b00 & 0b00000 = \phantom{2}0\cr
9 & 0b10 & 0b01 & 0b00100 = \phantom{2}4\cr
10 & 0b10 & 0b10 & 0b01000 = \phantom{2}8\cr
11 & 0b10 & 0b11 & 0b01100 = 12 \cr
12 & 0b11 & 0b00 & 0b00000 = \phantom{2}0\cr
13 & 0b11 & 0b01 & 0b01000 = \phantom{2}8\cr
14 & 0b11 & 0b10 & 0b10000 = 16 \cr
15 & 0b11 & 0b11 & 0b11000 = 24 \cr
\hline
\end{tabular}
\end{center}
It is clear this is neither monotonic nor efficient, as the same values appear
for several representations.
A better representation is made by realising that for all but the lowest
exponent representations, there is always a leading bit. Hence, this does not
have to be stored explicitly. This means this leading bit must be added to the
mantissa before the bit shift, and since this increments the length by one bit,
then the shift up needed is only $E-1$. The table below shows this improved
representation. It is monotonic, there are no duplicates, and the lowest two
exponent ranges give an exact representation.
\bigskip
\begin{center}
\begin{tabular}{c|c|c|c}
\hline
Representation & Exponent & Mantissa & Value \cr\hline
0 & 0b00 & 0b00 & 0b00000 = \phantom{2}0\cr
1 & 0b00 & 0b01 & 0b00001 = \phantom{2}1\cr
2 & 0b00 & 0b10 & 0b00010 = \phantom{2}2\cr
3 & 0b00 & 0b11 & 0b00011 = \phantom{2}3\cr
4 & 0b01 & 0b00 & 0b00100 = \phantom{2}4\cr
5 & 0b01 & 0b01 & 0b00101 = \phantom{2}5\cr
6 & 0b01 & 0b10 & 0b00110 = \phantom{2}6\cr
7 & 0b01 & 0b11 & 0b00111 = \phantom{2}7\cr
8 & 0b10 & 0b00 & 0b01000 = \phantom{2}8\cr
9 & 0b10 & 0b01 & 0b01010 = 10 \cr
10 & 0b10 & 0b10 & 0b01100 = 12 \cr
11 & 0b10 & 0b11 & 0b01110 = 14 \cr
12 & 0b11 & 0b00 & 0b10000 = 16 \cr
13 & 0b11 & 0b01 & 0b10100 = 20 \cr
14 & 0b11 & 0b10 & 0b11000 = 24 \cr
15 & 0b11 & 0b11 & 0b11100 = 28 \cr
\hline
\end{tabular}
\end{center}
In this improved representation, the mantissa has $M+1$ bits (except in
the lowest exponent range). The exponent can represent numbers up to
$2^E-1$ and hence will bit shift by a maximum of $2^E-2$ bits.
Hence, the number of bits in the representation is $M+E$ bits, while
the maximum number represented has $M+1+2^E-2 = M+2^E-1$ bits, i.e. is
is less than $2^{M+2^E-1}$. The reduction is $2^E-E-1$ bits.
For the example of $M=2$, $E=2$
in the table above, this gives $2-1+4=5$ bits, i.e.
numbers up to $2^5=32$ as shown and the reduction is 1 bit.
For the extreme values of $E$, then $E=0$ and $E=1$ both give an exact
representation as they only use the lowest range or two lowest ranges,
respectively. The reduction is $2^0-0-1=0$ and $2^1-1-1=0$ bits in both cases.
Explicitly, for $E=0$, then the representation has $M$ bits
while the value represented has $M$
bits also, i.e. the reduction is 0 bits.
For $E=1$, the representation has $M+1$ bits, while the
value represented has $M+1$ bits also, again with a reduction of 0 bits.
For the extreme value of $M=0$, then the representation is
just the number of bits in the input word. E.g. for $M=0$, $E=4$, then the
table is given below. The value range is less thn $2^{15}=32768$.
\bigskip
\begin{center}
\begin{tabular}{c|c|c|c}
\hline
Representation & Exponent & Mantissa & Value \cr\hline
0 & 0b0000 & 0 & 0b000000000000000 = \phantom{1222}0\cr
1 & 0b0001 & 0 & 0b000000000000001 = \phantom{1222}1\cr
2 & 0b0010 & 0 & 0b000000000000010 = \phantom{1222}2\cr
3 & 0b0011 & 0 & 0b000000000000100 = \phantom{1222}4\cr
4 & 0b0100 & 0 & 0b000000000001000 = \phantom{1222}8\cr
5 & 0b0101 & 0 & 0b000000000010000 = \phantom{122}16\cr
6 & 0b0110 & 0 & 0b000000000100000 = \phantom{122}32\cr
7 & 0b0111 & 0 & 0b000000001000000 = \phantom{122}64\cr
8 & 0b1000 & 0 & 0b000000010000000 = \phantom{12}128\cr
9 & 0b1001 & 0 & 0b000000100000000 = \phantom{12}256 \cr
10 & 0b1010 & 0 & 0b000001000000000 = \phantom{12}512 \cr
11 & 0b1011 & 0 & 0b000010000000000 = \phantom{1}1024 \cr
12 & 0b1100 & 0 & 0b000100000000000 = \phantom{1}2048 \cr
13 & 0b1101 & 0 & 0b001000000000000 = \phantom{1}4096 \cr
14 & 0b1110 & 0 & 0b010000000000000 = \phantom{1}8192 \cr
15 & 0b1111 & 0 & 0b100000000000000 = 16384 \cr
\hline
\end{tabular}
\end{center}
The maximum bit lengths of the value, i.e. $M+2^E-1$, for various values of
$M$ and $E$ are shown in the table below.
\bigskip
\begin{center}
\begin{tabular}{r||c|c|c|c|c|c|c|c|c}
\hline
$M=$ & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \cr\hline
$E=0$ & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \cr
1 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \cr
2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 \cr
3 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 \cr
4 & 15 & 16 & 17 & 18 & 19 & 20 & 21 & 22 & 23 \cr
5 & 31 & 32 & 33 & 34 & 35 & 36 & 37 & 38 & 39 \cr
6 & 63 & 64 & 65 & 66 & 67 & 68 & 69 & 70 & 71 \cr
7 & 127 & 128 & 129 & 130 & 131 & 132 & 133 & 134 & 135 \cr
8 & 255 & 256 & 257 & 258 & 259 & 260 & 261 & 262 & 263 \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
......
interface/AdcReading.hh
View file @ 4a65be59
......@@ -6,6 +6,8 @@
#include <cassert>
#include <cmath>
#include "Random.hh"
class AdcReading {
public:
AdcReading() : fLoData(0), fHiData(0) {
......@@ -14,6 +16,11 @@ public:
~AdcReading() {
}
void initialise() {
fLoGain=Random::random().Gaus();
fHiGain=Random::random().Uniform()-0.5;
}
uint16_t loData() const {
return fLoData;
}
......@@ -33,13 +40,15 @@ public:
fHiData=0;
} else {
double loAdcLsb(fLoAdcLsb*(1.0+fLoGainWidth*fLoGain));
double hiAdcLsb(fHiAdcLsb*(1.0+fHiGainWidth*fHiGain));
// TP = baseline
if(fArchitecture==0) {
// Low range
fLoData=1023;
if(q<1023.0*fLoAdcLsb) fLoData=uint16_t(q/fLoAdcLsb);
if(q<1023.0*loAdcLsb) fLoData=uint16_t(q/loAdcLsb);
// High range
if(q>100.0) q-=100.0*(100.0-90.0)/(q-90.0);
......@@ -54,7 +63,7 @@ public:
// Low range
fLoData=1023;
if(q<1023.0*fLoAdcLsb) fLoData=uint16_t(q/fLoAdcLsb);
if(q<1023.0*loAdcLsb) fLoData=uint16_t(q/loAdcLsb);
// High range
fHiData=16383;
......@@ -68,12 +77,42 @@ public:
// Low range
fLoData=255;
if(q< 255.0*fLoAdcLsb) fLoData=uint16_t(q/fLoAdcLsb);
if(q< 255.0*loAdcLsb) fLoData=uint16_t(q/loAdcLsb);
// High range
fHiData=4095;
if(q<4095.0*fHiAdcLsb) fHiData=uint16_t(q/fHiAdcLsb);
}
// Single gain ~ideal
if(fArchitecture==3) {
// Low and high ranges combined to store more than 16 bits
uint32_t total(524287);
if(q<524287.0*loAdcLsb) total=uint32_t(q/loAdcLsb);
fLoData=(total&0xffff);
fHiData=(total>>16);
if(fHiData>7) {
std::cout << "AdcReading::reading() Arch 3, q = " << q
<< ", loAdcLsb = " << loAdcLsb
<< ", fLoData = " << fLoData
<< ", fHiData = " << fHiData
<< ", total = " << total << std::endl;
}
}
// Bi-gain ~ideal
if(fArchitecture==4) {
// Low range
fLoData=1023;
if(q<1023.0*loAdcLsb) fLoData=uint16_t(q/loAdcLsb);
// High range
fHiData=4095;
if(q<4095.0*hiAdcLsb) fHiData=uint16_t(q/hiAdcLsb);
}
}
}
......@@ -143,6 +182,22 @@ public:
return convertLsb(fHiData,fHiAdcLsb);
}
if(fArchitecture==3) {
uint32_t total(fHiData);
total<<=16;
/*
std::cout << "AdcReading::result() Arch 3, fLoData = " << fLoData
<< ", fHiData = " << fHiData
<< ", result = " << total+fLoData << std::endl;
*/
return total+fLoData;
}
if(fArchitecture==4) {
if(fLoData<1023) return fLoData;
return 25*fHiData+12;
}
assert(false);
return 0xffffffff;
}
......@@ -163,7 +218,12 @@ public:
assert(false);
return 999999.0;
*/
return fAdcLsb*(result()+0.5);
unsigned r(result());
if(fTruncateBits>0) {
r>>=fTruncateBits;
r<<=fTruncateBits;
}
return fAdcLsb*(r+0.5);
}
bool read(std::istream &fin) {
......@@ -196,6 +256,8 @@ public:
if(fArchitecture==0) return 5;
if(fArchitecture==1) return 8;
if(fArchitecture==2) return 5;
if(fArchitecture==3) return 1;
if(fArchitecture==4) return 1;
return 0;
}
......@@ -230,6 +292,20 @@ public:
fRatio=50;
fLoAdcLsb=fHiAdcLsb/fRatio;
}
if(fArchitecture==3) {
fMethod=0;
fLoAdcLsb=fAdcLsb;
fRatio=0.0;
fHiAdcLsb=0.0;
}
if(fArchitecture==4) {
fMethod=0;
fLoAdcLsb=fAdcLsb;
fRatio=25.0;
fHiAdcLsb=fRatio*fAdcLsb;
}
}
static const unsigned fNumberOfArchitectures;
......@@ -241,15 +317,22 @@ public:
static double fLoAdcLsb;
static double fHiAdcLsb;
static double fLoGainWidth;
static double fHiGainWidth;
static unsigned fTruncateBits;
protected:
static unsigned fArchitecture;
static unsigned fMethod;
uint16_t fLoData;
uint16_t fHiData;
float fLoGain;
float fHiGain;
};
const unsigned AdcReading::fNumberOfArchitectures=3;
const unsigned AdcReading::fNumberOfArchitectures=5;
unsigned AdcReading::fArchitecture=0;
const unsigned AdcReading::fMaximumNumberOfMethods=8;
......@@ -261,4 +344,9 @@ double AdcReading::fLoAdcLsb=AdcReading::fAdcLsb;
double AdcReading::fRatio=24.9;
double AdcReading::fHiAdcLsb=AdcReading::fRatio*AdcReading::fAdcLsb;
double AdcReading::fLoGainWidth=0.01;
double AdcReading::fHiGainWidth=1/25.0;
unsigned AdcReading::fTruncateBits=0;
#endif
interface/AnalysisPhoton.hh
View file @ 4a65be59
......@@ -19,6 +19,7 @@
#include "AnalysisBase.hh"
#include "SubAnalysisPhotonSimHit.hh"
#include "SubAnalysisPhotonTrgHit.hh"
#include "SubAnalysisPhotonFeArchitecture.hh"
#include "Event.hh"
#include "ErrorMatrix.hh"
......@@ -244,13 +245,20 @@ public:
AnalysisPhoton(const std::string &sRoot="", const std::string &sOut="")
: AnalysisBase("AnalysisPhoton",sRoot,sOut),
fSAPSimHit(sRoot,sOut), fSAPTrgHit(sRoot,sOut) {
fSAPSimHit(sRoot,sOut),
fSAPTrgHit(sRoot,sOut),
fSAPFeArchitecture(sRoot,sOut) {
// What to do?
fSimEvent=true;
fSimEvent=false;
fDigEvent=false;
fTrgEvent=true;
//fTrgEvent=true;
fTrgEvent=false;
fFeArch=true;
fSelect=false;
fC2DEvent=false;
fC3DEvent=false;
......@@ -2732,6 +2740,7 @@ public:
}
virtual bool select(Event &event) {
if(!fSelect) return true;
bool good(false);
for(unsigned e(0);e<Geometry::kNumberOfEndcaps;e++) {
......@@ -2846,6 +2855,7 @@ public:
//event->trgHitProcessor();
if(fTrgEvent) fSAPTrgHit.event(event);
if(fFeArch) fSAPFeArchitecture.event(event);
return true;
for(unsigned e(0);e<Geometry::kNumberOfEndcaps;e++) {
......@@ -3371,6 +3381,8 @@ protected:
bool fSimEvent;
bool fDigEvent;
bool fTrgEvent;
bool fFeArch;
bool fSelect;
bool fC2DEvent;
bool fC3DEvent;
......@@ -3378,6 +3390,7 @@ protected:
SubAnalysisPhotonSimHit fSAPSimHit;
SubAnalysisPhotonTrgHit fSAPTrgHit;
SubAnalysisPhotonFeArchitecture fSAPFeArchitecture;
TH1F *hSelectSinThetaAll,*hSelectSinThInvAll;
......
interface/Coefficients.hh
View file @ 4a65be59
......@@ -55,7 +55,7 @@ public:
}
double solve(uint64_t lm=0xffffffffffffffff) {
assert(fNum>0.5);
if(fNum<0.5) return -999.0;
assert(lm!=0);
uint64_t one(1);
......@@ -148,6 +148,11 @@ public:
return (fErg/fNum)-sum*fCoe;
}
double energy(const TVectorD &p) const {
assert(p.GetNrows()==kNumberOfCoefficients);
return fCoe*p;
}
const TMatrixDSym& coefficientMatrix() const {
return fSSq;
}
......
interface/DigHit.hh
View file @ 4a65be59
......@@ -80,19 +80,25 @@ public:
vSimHit=&s;
return true;
}
double sigma() const {
return DigNoise::noise(Geometry::point(*this).rho(),
Geometry::numberOfSublayers(layer(),wafer())-1,
Geometry::halfCell(*this));
}
void setNoise() {
fGaussian=Random::random().Gaus();
}
double noise() const {
return fGaussian*DigNoise::noise(Geometry::point(*this).rho(),
Geometry::numberOfSublayers(layer(),wafer())-1,
Geometry::halfCell(*this));
return fGaussian*sigma();
}
void process() {
reading(simCharge()+noise());
reading(inputCharge());
}
// Result
......@@ -100,12 +106,20 @@ public:
return chargeResult();
}
double sigmas() const {
return charge()/sigma();
}
double deposit() const {
return chargeResult()/Constants::fFcPerGev;
return charge()/Constants::fFcPerGev;
}
double mips() const {
return chargeResult()*Constants::fMipPerFc[Geometry::numberOfSublayers(layer(),wafer())-1];
return charge()*Constants::fMipPerFc[Geometry::numberOfSublayers(layer(),wafer())-1];
}
double transverseMips() const {
return mips()*Geometry::point(*this).sinTheta();
}
double energy() const {
......@@ -116,6 +130,24 @@ public:
return energy()*Geometry::point(*this).sinTheta();
}
// Access to input charge (pre-digitisation) information
double inputCharge() const {
return simCharge()+noise();
}
double inputSigmas() const {
return inputCharge()/sigma();
}
double inputMips() const {
return inputCharge()*Constants::fMipPerFc[Geometry::numberOfSublayers(layer(),wafer())-1];
}
double inputTransverseMips() const {
return inputMips()*Geometry::point(*this).sinTheta();
}
// Access to SimHit information
unsigned numberOfSimHits() const {
......
interface/Event.hh
View file @ 4a65be59
......@@ -135,6 +135,20 @@ public:
return q;
}
double simTransverseMips(unsigned e, unsigned l, unsigned w, unsigned c, bool signal=false) const {
assert(e<Geometry::kNumberOfEndcaps);
assert(l<Geometry::kNumberOfHGCLayers);
assert(Geometry::validWafer(l,w));
assert(c<Geometry::numberOfCells(l,w));
double q(0.0);
for(unsigned i(0);i<vSimHit[e][l][w][c].size();i++) {
if(!signal || (signal && vSimHit[e][l][w][c][i].isSignal()))
q+=vSimHit[e][l][w][c][i].transverseMips();
}
return q;
}
double simTransverseEnergy(unsigned e, unsigned l, unsigned w, unsigned c, bool signal=false) const {
assert(e<Geometry::kNumberOfEndcaps);
assert(l<Geometry::kNumberOfHGCLayers);
......@@ -641,6 +655,7 @@ public:
if(Geometry::validWafer(l,w)) {
for(unsigned c(0);c<Geometry::numberOfCells(l,w);c++) {
vDigHit[e][l][w][c].process();
vDigHit[e][l][w][c].selected(true); // SHOULD ADD CUT IN SIGMAS!
}
}
}
......@@ -701,6 +716,7 @@ public:
for(unsigned c(0);c<Geometry::numberOfCells(l,w);c++) {
vDigHit[e][l][w][c].detId(e,l,w,Geometry::waferType(l,w),c);
vDigHit[e][l][w][c].initialise();
vDigHit[e][l][w][c].addSimHitVector(vSimHit[e][l][w][c]);
vTrgHit[e][l][w][Geometry::triggerCell(vDigHit[e][l][w][c])].addDigHit(vDigHit[e][l][w][c]);
}
......
interface/Geometry.hh
View file @ 4a65be59
......@@ -322,7 +322,7 @@ public:
std::cout << std::setw(4) << vBcToTc128[i][j];
Point2D xy(cellPoint2D(true,vBcToTc128[i][j]));
//std::cout << " " << xy.x() << " " << xy.y();
std::cout << " " << xy.x() << " " << xy.y();
w=(halfCell(true,vBcToTc128[i][j])?0.5:1.0);
wSum+=w;
......@@ -1373,6 +1373,38 @@ public:
}
}
static uint32_t cellFlip(bool t128, uint32_t c) {
std::pair<int,int> xy(cellPair(t128,c));
if(xy.first==0) {
//std::cout << "Geometry::cellFlip() 128-cell sensor cell " << c
// << " not flipped" << std::endl;
return c;
}
xy.first=-xy.first;
if(t128) {
for(unsigned i(0);i<kNumberOfType0Cells;i++) {
if(cellPair(t128,i)==xy) {
//std::cout << "Geometry::cellFlip() 128-cell sensor cell " << c
// << " flipped to cell " << i << std::endl;
return i;
}
}
} else {
for(unsigned i(0);i<kNumberOfType1Cells;i++) {
if(cellPair(t128,i)==xy) {
//std::cout << "Geometry::cellFlip() 256-cell sensor cell " << c
// << " flipped to cell " << i << std::endl;
return i;
}
}
}
assert(false);
return 0xffffffff;
}
/////////////////////////////////////////////////////////////////
// Dimensions in 2D
......@@ -1479,7 +1511,7 @@ public:
<< cellPoint2D(waferType(l,w)==0,c).y() << std::endl;
*/
Point2D pxy(cellPoint2D(waferType(l,w)==0,c));
if(e==0) pxy.reflectX();
//if(e==0) pxy.reflectX(); Fixed in ReadSimHitBase
return Point(waferPoint(e,l,w)+pxy);
}
......