Commit 58a3450b authored by dauncey's avatar dauncey

Add info on pseudorapidity

parent 392a2e92
void Rapidity() {
gROOT->SetStyle("Plain");
TCanvas *c1(new TCanvas());
gStyle->SetOptStat(0);
gStyle->SetPalette(1);
//gStyle->SetPalette(kDarkBodyRadiator);
//gStyle->SetPalette(kBlueGreenYellow);
/*
cout << c1->GetLeftMargin();
c1->SetLeftMargin(0.1);
cout << " to " << c1->GetLeftMargin() << endl;
cout << c1->GetRightMargin();
c1->SetRightMargin(0.14);
cout << " to " << c1->GetRightMargin() << endl;
*/
TLatex t;
TH1F *h1;
TH2F *h2[10];
h2[0]=new TH2F("H20",";#eta;sin #theta",16,1.4,3.0,1,0.0,0.5);
TF1 st("st","sin(atan(1/sinh(x)))",1.4,3.0);
h2[0]->Draw();
st.Draw("same");
c1->Print("SinThetaVsEta.pdf");
h2[1]=new TH2F("H21",";#eta;Derivative",16,1.4,3.0,1,0.0,10.0);
TF1 deta("deta","tanh(x)*sinh(x)",1.4,3.0);
TF1 dphi("dphi","sinh(x)",1.4,3.0);
h2[1]->Draw();
deta.Draw("same");
dphi.SetLineColor(2);
dphi.Draw("same");
t.SetTextColor(1);
t.DrawLatex(1.5,9.0,"#partial#eta/#partial#rho_{s}");
t.SetTextColor(2);
t.DrawLatex(1.5,8.0,"#partial#phi/#partials_{s}");
c1->Print("DerivativesVsEtaA.pdf");
h2[2]=new TH2F("H22",";#eta;Derivative",16,1.4,3.0,6,0.0,0.6);
TF1 drho("drho","1/(tanh(x)*sinh(x))",1.4,3.0);
TF1 ds("ds","1/sinh(x)",1.4,3.0);
h2[2]->Draw();
drho.Draw("same");
ds.SetLineColor(2);
ds.Draw("same");
t.SetTextColor(1);
t.DrawLatex(2.7,0.50,"#partial#rho_{s}/#partial#eta");
t.SetTextColor(2);
t.DrawLatex(2.7,0.45,"#partials_{s}/#partial#phi");
c1->Print("DerivativesVsEtaB.pdf");
h2[3]=new TH2F("H23",";#eta;1/Jacobian",16,1.4,3.0,1,0.0,100.0);
TF1 dj("dj","tanh(x)*sinh(x)*sinh(x)",1.4,3.0);
h2[3]->Draw();
dj.Draw("same");
c1->Print("JacobianVsEta.pdf");
}
No preview for this file type
\documentclass[10pt]{article}
\usepackage{a4}
\usepackage{amsmath}
\usepackage{graphicx}
\oddsidemargin=0pt % No extra space wasted after first inch.
\evensidemargin=0pt % Ditto (for two-sided output).
......@@ -216,6 +217,15 @@ around the truth value is given by
\mathrm{RMS}^2 = \frac{1}{N} \sum_e (E_e - T_e)^2
= \frac{1}{N} \sum_e \left(\sum_i a_i d_{e,i} - T_e\right)^2
\end{equation}
Alternatively, we could optimise for $\Delta E/E$, i.e.
\begin{equation}
\mathrm{RMS}^2 = \frac{1}{N} \sum_e \left(\frac{E_e - T_e}{T_e}\right)^2
= \frac{1}{N} \sum_e \left(\frac{\sum_i a_i d_{e,i} - T_e}{T_e}\right)^2
= \frac{1}{N} \sum_e \left(\sum_i a_i d_{e,i}^\prime - 1\right)^2
\end{equation}
where $d_{e,i}^\prime = d_{e,i}/T_e$. Hence, this method will be identical in form
in the following, but replacing $T_e$ by 1 and $d_{e,i}$ by $d_{e,i}^\prime$.
This can be thought of as similar to a chi-squared; we want to minimise this
expression. If all the $a_i$ are considered as independent parameters (so 28 for
the EE only), then explicitly
......@@ -466,6 +476,29 @@ 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.
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}
%\newpage
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.
......@@ -504,29 +537,6 @@ Representation & Exponent & Mantissa & Value \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}
\newpage
\section{Template fit of energy in depth}
Assume a template shape for a photon of $P_l$ per photon energy GeV
for layer $l$.
......@@ -708,6 +718,138 @@ Again, for small $\Delta\phi$, then $\overline{\rho}$
clearly approximates to the geometric
mean of the two radii, as before.
\section{Cartesian coordinates and pseudorapidity}
As above, take $\rho = \sqrt{x^2+y^2}$ and
define the arc length $s=\rho\phi$. In terms of scaled coordinates $x_s=x/z$
and $y_s=y/z$, then equivalently $\rho_s = \rho/z$ and
$s_s = (\rho/z)\phi = \rho_s\phi$. Note, $\rho_s = \tan\theta$ and so
$x_s = \tan\theta\cos\phi$ and $y_s = \tan\theta\sin\phi$.
Generally $\sinh\eta = z/\rho = 1/\rho_s$, so $\rho_s = 1/\sinh\eta$
and $s_s = \phi/\sinh\eta$. To convert from $E$ to $E_T$ requires
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]
\end{equation}
%
\begin{figure}[ht!]
\begin{center}
\includegraphics[width=9cm]{SinThetaVsEta.pdf}
%\caption{
% Left: Active rotation by $30^\circ$
% of an object in a fixed coordinate system.
% Right: Passive rotation by $30^\circ$
% of the coordinate system around a fixed object.
%}
\label{fig:rotations}
\end{center}
\end{figure}
Comparing areas in $\Delta \rho_s \Delta s_s$ and $\Delta \eta \Delta \phi$
requires the Jacobian. From above
\begin{equation}
\frac{\partial \rho_s}{\partial \eta} = - \frac{\cosh\eta}{\sinh^2\eta},\qquad
\frac{\partial \rho_s}{\partial \phi} = 0
\end{equation}
while
\begin{equation}
\frac{\partial s_s}{\partial \eta} = - \frac{\phi\cosh\eta}{\sinh^2\eta},\qquad
\frac{\partial s_s}{\partial \phi} = \frac{1}{\sinh\eta}
\end{equation}
Hence the Jacobian is
\begin{equation}
J = \left|\frac{\partial \rho_s}{\partial \eta}\frac{\partial s_s}{\partial \phi}
- \frac{\partial \rho_s}{\partial \phi}\frac{\partial s_s}{\partial \eta}\right|
= \left|\frac{\cosh\eta}{\sinh^2\eta}\times\frac{1}{\sinh\eta}\right|
= \left|\frac{\cosh\eta}{\sinh^3\eta}\right|
\end{equation}
Hence
\begin{equation}
\Delta \rho_s \Delta s_s = \left|\frac{\cosh\eta}{\sinh^3\eta}\right|
\Delta \eta \Delta \phi
\end{equation}
Cross-check the other way around, using
\begin{equation}
\eta = \sinh^{-1}\left(\frac{1}{\rho_s}\right),\qquad
\phi = \frac{s_s}{\rho_s}
\end{equation}
then
\begin{equation}
\frac{\partial \eta}{\partial \rho_s}
= - \frac{\sinh^2\eta}{\cosh\eta} = -\frac{1}{\rho_s^2\sqrt{1+1/\rho_s^2}}
= -\frac{1}{\rho_s\sqrt{1+\rho_s^2}},\qquad
\frac{\partial \eta}{\partial s_s} = 0
\end{equation}
and
\begin{equation}
\frac{\partial \phi}{\partial \rho_s}
= - \frac{s_s}{\rho_s^2},\qquad
\frac{\partial \phi}{\partial s_s} = \frac{1}{\rho_s}
\end{equation}
Hence the inverse Jacobian to the above is
\begin{equation}
\frac{1}{J} = \left|\frac{\partial \eta}{\partial \rho_s}\frac{\partial \phi}{\partial s_s}
-\frac{\partial \phi}{\partial \rho_s}\frac{\partial \eta}{\partial s_s}\right|
= \left|\frac{1}{\rho_s\sqrt{1+\rho_s^2}}\times\frac{1}{\rho_s}\right|
= \left|\frac{1}{\rho_s^2\sqrt{1+\rho_s^2}}\right|
\end{equation}
In terms of $\eta$, this is
\begin{equation}
\frac{1}{J}
= \left|\frac{1/\rho_s^3}{\sqrt{1+1/\rho_s^2}}\right|
= \left|\frac{\sinh^3\eta}{\cosh\eta}\right|
\end{equation}
and so is consistent. The derivatives are shown below, showing that they
are very similar in size, so that a circle in $\rho_s,s_s$ coordinates
is also circular to a good approximation in $\eta,\phi$. At the innermost
edge of the HGCAL, i.e. $\eta=3$, the derivatives are around 10, so for
example $\Delta\rho_s=\Delta s_s=0.02$ is approximately equivalent to
$\Delta\eta=\Delta\phi = 0.2$.
%
\begin{figure}[ht!]
\begin{center}
\includegraphics[width=7cm]{DerivativesVsEtaA.pdf}
\hspace{1cm}
\includegraphics[width=7cm]{DerivativesVsEtaB.pdf}
%\caption{
% Left: Active rotation by $30^\circ$
% of an object in a fixed coordinate system.
% Right: Passive rotation by $30^\circ$
% of the coordinate system around a fixed object.
%}
\label{fig:derivatives}
\end{center}
\end{figure}
The size of the inverse Jacobian is shown below, which gives the factor for
converting areas. Specifically, with a pileup transverse energy density
at PU 200 of
around 140\,GeV per unit area in $\Delta\eta\Delta\phi$, then the average
transverse energy in GeV in a circle of $\Delta\rho_s=\Delta s_s=0.02$ is
approximately equivalent to
\begin{equation}
E_T = \frac{140}{J} 0.02^2 \pi
\end{equation}
At $\eta=1.4$, this is around 0.5\,GeV, while
at $\eta=3.0$, this is around 18\,GeV.
%
\begin{figure}[ht!]
\begin{center}
\includegraphics[width=9cm]{JacobianVsEta.pdf}
%\caption{
% Left: Active rotation by $30^\circ$
% of an object in a fixed coordinate system.
% Right: Passive rotation by $30^\circ$
% of the coordinate system around a fixed object.
%}
\label{fig:jacobian}
\end{center}
\end{figure}
\section{Shower position and direction fit}
\section{Motion in a magnetic field}
......
......@@ -9,18 +9,39 @@
class Coefficients {
public:
enum {
kNumberOfCoefficients=Geometry::kNumberOfLayers
kNumberOfLayers=Geometry::kNumberOfLayers,
kNumberOfCoefficients=kNumberOfLayers+1
};
Coefficients(bool b=false) : fPlainRms(b),
Coefficients(unsigned n=0, bool c=false) :
fSum(kNumberOfCoefficients),
fXSum(kNumberOfCoefficients),
fCoe(kNumberOfCoefficients),
fSSq(kNumberOfCoefficients) {
initialise(n,c);
}
Coefficients(bool b, bool c=false) :
fSum(kNumberOfCoefficients),
fXSum(kNumberOfCoefficients),
fCoe(kNumberOfCoefficients),
fSSq(kNumberOfCoefficients) {
initialise(b?0:1,c);
}
void initialise(unsigned n, bool c) {
std::cout << std::endl << "Coefficients::initialise() called with power = "
<< n << ", constant = " << (c?"true":"false") << std::endl;
fPlainRms=n;
fConstant=c;
fNum=0.0;
fErg=0.0;
fXErg=0.0;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fSum(i)=0.0;
fXSum(i)=0.0;
fCoe(i)=0.0;
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
fSSq(i,j)=0.0;
......@@ -30,30 +51,47 @@ public:
}
void fill(double e, const double *p) {
TVectorD v(kNumberOfCoefficients);
for(unsigned i(0);i<kNumberOfCoefficients;i++) v(i)=p[i];
TVectorD v(kNumberOfLayers);
for(unsigned i(0);i<kNumberOfLayers;i++) v(i)=p[i];
fill(e,v);
}
void fill(double e, const TVectorD &p) {
assert(p.GetNrows()==kNumberOfCoefficients);
assert(p.GetNrows()==kNumberOfLayers);
fNum++;
if(fPlainRms) fErg+=e*e;
else fErg+=e;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
if(fPlainRms) fSum(i)+=e*p(i);
else fSum(i)+=p(i);
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
if(fPlainRms) fSSq(i,j)+=p(i)*p(j);
else fSSq(i,j)+=p(i)*p(j)/e;
double k(1);
if(fPlainRms==1) k=1.0/(e);
if(fPlainRms==2) k=1.0/(e*e);
if(fPlainRms==3) k=1.0/(e*e*e);
if(fPlainRms==4) k=1.0/(e*e*e*e);
if(fPlainRms==5) k=1.0/(e*e*e*e*e);
if(fPlainRms==6) k=1.0/(e*e*e*e*e*e);
fErg+=e*e*k;
fXErg+=e*sqrt(k);
for(unsigned i(0);i<kNumberOfLayers;i++) {
fSum(i)+=e*p(i)*k;
fXSum(i)+=p(i)*sqrt(k);
for(unsigned j(0);j<kNumberOfLayers;j++) {
fSSq(i,j)+=p(i)*p(j)*k;
}
}
fSum(kNumberOfLayers)+=e*k;
fXSum(kNumberOfLayers)+=sqrt(k);
for(unsigned j(0);j<kNumberOfLayers;j++) {
fSSq(kNumberOfLayers,j)+=p(j)*k;
fSSq(j,kNumberOfLayers)+=p(j)*k;
}
fSSq(kNumberOfLayers,kNumberOfLayers)+=k;
}
void fillCoe(const double *c) {
for(unsigned i(0);i<kNumberOfCoefficients;i++) fCoe(i)=c[i];
}
double solve(uint64_t lm=0xffffffffffffffff) {
if(fNum<0.5) return -999.0;
assert(lm!=0);
......@@ -61,14 +99,21 @@ public:
uint64_t one(1);
TVectorD sum(fSum);
TVectorD xsum(fXSum);
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fCoe(i)=0.0;
bool ignore((lm&(one<<i))==0);
if(ignore) sum(i)=0.0;
else sum(i)/=fNum;
else {
sum(i)/=fNum;
xsum(i)/=fNum;
}
}
if(!fConstant) sum(kNumberOfCoefficients-1)=0.0;
TMatrixDSym ssq(fSSq);
//TMatrixD ssq(fSSq);
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
......@@ -143,14 +188,23 @@ public:
<< inv(0,1) << std::endl;
std::cout << "solve fErg = " << fErg/fNum << ", sum*fCoe = "
<< sum*fCoe << std::endl;
<< sum*fCoe << ", chi-sq = "
<< (fErg/fNum)-sum*fCoe << std::endl;
std::cout << "solve fXErg = " << fXErg/fNum << ", xsum*fCoe = "
<< xsum*fCoe << ", mean = "
<< (fXErg/fNum)-xsum*fCoe << std::endl;
return (fErg/fNum)-sum*fCoe;
}
double energy(const TVectorD &p) const {
assert(p.GetNrows()==kNumberOfCoefficients);
return fCoe*p;
assert(p.GetNrows()==kNumberOfLayers);
double e(fCoe(kNumberOfLayers));
for(unsigned i(0);i<kNumberOfLayers;i++) e+=fCoe(i)*p(i);
return e;
}
const TMatrixDSym& coefficientMatrix() const {
......@@ -161,6 +215,16 @@ public:
return fCoe;
}
void operator+=(const Coefficients &c) {
assert(fPlainRms==c.fPlainRms);
fNum+=c.fNum;
fErg+=c.fErg;
fSum+=c.fSum;
fSSq+=c.fSSq;
for(unsigned i(0);i<kNumberOfCoefficients;i++) fCoe(i)=0.0;
}
bool read(const std::string &fName) {
std::ifstream fin(fName.c_str());
if(!fin) return false;
......@@ -188,16 +252,6 @@ public:
return true;
}
void operator+=(const Coefficients &c) {
assert(fPlainRms==c.fPlainRms);
fNum+=c.fNum;
fErg+=c.fErg;
fSum+=c.fSum;
fSSq+=c.fSSq;
for(unsigned i(0);i<kNumberOfCoefficients;i++) fCoe(i)=0.0;
}
bool write(const std::string &fName) {
std::ofstream fout(fName.c_str());
if(!fout) return false;
......@@ -225,6 +279,42 @@ public:
return true;
}
bool readCoe(const std::string &fName) {
std::ifstream fin(fName.c_str());
if(!fin) return false;
unsigned n;
fin >> n;
if(n!=kNumberOfCoefficients) return false;
Double_t d;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fin >> d;
fCoe(i)=d;
}
if(!fin) return false;
return true;
}
bool writeCoe(const std::string &fName) {
std::ofstream fout(fName.c_str());
if(!fout) return false;
fout << kNumberOfCoefficients
<< std::endl;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fout << " " << std::setprecision(20) << fCoe(i);
}
fout << std::endl;
if(!fout) return false;
fout.close();
return true;
}
void print() {
std::cout << "Coefficients::print() Number of entries = "
<< fNum << ", energy^2 " << fErg << std::endl;
......@@ -249,10 +339,15 @@ public:
std::cout << std::endl;
}
const bool fPlainRms;
unsigned fPlainRms;
bool fConstant;
double fNum,fErg;
TVectorD fSum,fCoe;
TMatrixDSym fSSq;
TVectorD fXSum;
double fXErg;
};
#endif
......@@ -11,6 +11,10 @@
#include "TH1F.h"
#include "TLorentzVector.h"
#include "HepMC/GenEvent.h"
#include "Point.hh"
class TrgC3d;
class Particle {
......@@ -91,7 +95,13 @@ public:
fFourP=TLorentzVector( p4.x(), p4.y(), p4.z(), p4.e() );
}
Particle(int pdg, double px, double py, double pz) : pGenParticle(0) {
fPdgId=pdg;
fPosition=Point(0,0,0);
fMomentum=Point(px,py,pz);
}
~Particle() {
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment