diff --git a/doc/ExponentialCumulative.pdf b/doc/ExponentialCumulative.pdf new file mode 100644 index 0000000000000000000000000000000000000000..18b9eb39090b0ec5e88f45eb8a4e90f1c34f4a00 Binary files /dev/null and b/doc/ExponentialCumulative.pdf differ diff --git a/doc/ExponentialSqrtCumulative.pdf b/doc/ExponentialSqrtCumulative.pdf new file mode 100644 index 0000000000000000000000000000000000000000..b552103a10f67a6e62672dc1b5c0dccc1362d9aa Binary files /dev/null and b/doc/ExponentialSqrtCumulative.pdf differ diff --git a/doc/GaussianCumulative.pdf b/doc/GaussianCumulative.pdf new file mode 100644 index 0000000000000000000000000000000000000000..ff9add23646a9c492144b3776a8214cec96227c8 Binary files /dev/null and b/doc/GaussianCumulative.pdf differ diff --git a/doc/wup.pdf b/doc/wup.pdf index cb71a9454adb6c85d00d192e4a3494ff7b796faa..06ccd5efb7e6bdf38abb20e488a7227c4ee221f7 100644 Binary files a/doc/wup.pdf and b/doc/wup.pdf differ diff --git a/doc/wup.tex b/doc/wup.tex index d046410e0611cd63aaf0847ba0488529db742daa..efc71585086c281f6242ee5ed228204e4c7d55ef 100644 --- a/doc/wup.tex +++ b/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} diff --git a/interface/PackerBase.hh b/interface/PackerBase.hh index 37a72d5988a65843c2d1b28357701fd72222c3f4..5a7cc7310518d0297ecd62ba069a99892ceb6ece 100644 --- a/interface/PackerBase.hh +++ b/interface/PackerBase.hh @@ -8,24 +8,27 @@ #include "utilities.cc" -class PackerBase { +template +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 PackerBase; + #endif diff --git a/interface/PackerFloat22.hh b/interface/PackerFloat22.hh new file mode 100644 index 0000000000000000000000000000000000000000..7a1f1402be8050cc9247ef1393dfe622acdca9ef --- /dev/null +++ b/interface/PackerFloat22.hh @@ -0,0 +1,104 @@ +#ifndef PackerFloat22_HH +#define PackerFloat22_HH + +#include +#include +#include +#include + +#include "PackerBase.hh" + +class PackerFloat22 : public PackerBaseUint { +public: + PackerFloat22(unsigned e, unsigned m) : fMan(m), fExp(e), fBit(fMan+(1<>=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)<>fMan); + assert(e<(uint32_t(1)<>(16-e)); + } + + uint32_t dataMax() const { + return fDataMax; + } + + void print() const { + std::cout << "PackerFloat22::print()";PackerBaseUint::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 diff --git a/interface/TE1F.hh b/interface/TE1F.hh index 9ed840fc08ae335c764d7a1ca9ac97b23296a800..ad1398275aa7b9ab7011a2e92904fd065df2db1f 100644 --- a/interface/TE1F.hh +++ b/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);bFill(i->GetBinCenter(b),w); - } - if(x>i->GetBinCenter(n)) i->Fill(i->GetBinCenter(n),w); + if(up) { + for(int b(1);bFill(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(xGetBinCenter(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; diff --git a/src/DaqDaqBuffer.cpp b/src/DaqDaqBuffer.cpp index bb3cf28f41d8a1e5beebda5860c42c48aafabd25..af637b63775f202a7ad19cc281174afea955df71 100644 --- a/src/DaqDaqBuffer.cpp +++ b/src/DaqDaqBuffer.cpp @@ -1,3 +1,10 @@ +/* + + 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 diff --git a/src/Langford.cpp b/src/Langford.cpp new file mode 100644 index 0000000000000000000000000000000000000000..5591b98ef324486790f0634f1024565bf0592a0f --- /dev/null +++ b/src/Langford.cpp @@ -0,0 +1,741 @@ +#include "Hack.hh" + +#include +#include +#include +#include +#include + +#include "TFileHandler.hh" +#include "TH2B.hh" + +/* +NUMBERS OF HGCROCS PER MODULE FOR VARIOUS PARTIAL TYPES + +Full (type F): 6 ROCs for High Density, 3 for Low Density + +Five (type b): 5 ROCs for High Density, 3 for Low Density +Choptwo (type g) 5 ROCs for HD, 3 for LD +Halves (type a) 3 ROCs for HD, 2 for LD + +Then some guess +Semis (type d) 4 ROCs for HD, 2 for LD. -> CHANGE TO 3 FOR HD!!! +Thirds (type c) hopefully only one ROC for HD, 1 for LD + +The remaining ones should exist only in LD +Chopfour (h) 1 ROC +Type i (or y) hopefully only 1 ROC +Type j (or z) hopefully 2 ROCs +*/ + +class Hgcroc { +public: + Hgcroc() { + } + + void read(std::istream &fin) { + fin >> fMean >> fRms; + assert(fin); + print(); + } + + void print() const { + std::cout << "Hgcroc::print() Channels " << fChannels + << ", mean, RMS = " << fMean << ", " << fRms << std::endl; + } + + unsigned fChannels; + double fMean; + double fRms; +}; + +class Module { +public: + Module() { + } + + void read(std::istream &fin) { + char bra,ket,comma; + std::string sDensity; + + fin >> bra >> fU >> comma >> fV >> ket >> fShape >> fType >> sDensity; + assert(bra='('); + assert(ket=')'); + assert(comma=','); + type(); + numberOfHgcrocs(); + + unsigned nH(0); + if(sDensity=="high") nH=6; + if(sDensity=="low") nH=3; + assert(nH>0); + + unsigned mh; + for(unsigned h(0);h> mh; + assert(mh==h); + + vHgcroc.push_back(Hgcroc()); + vHgcroc[vHgcroc.size()-1].fChannels=(sDensity=="high"?432:192); + vHgcroc[vHgcroc.size()-1].read(fin); + } + assert(fin); + print(); + } + + unsigned type() const { + if(fType=='I') return 0; + if(fType=='M') return 1; + if(fType=='O') return 2; + std::cout << "Module::type() fType = " << fType << std::endl; + assert(false); + } + + unsigned numberOfDaqElinks() const { + return 2*numberOfHgcrocs(); + } + + unsigned numberOfTpgElinks() const { + if(fType=='I') return 2*numberOfHgcrocs(); + return 4*numberOfHgcrocs(); + } + + unsigned numberOfHgcrocs() const { + if(fType=='I') { + if(fShape=='F') return 6; + if(fShape=='b') return 5; + if(fShape=='g') return 5; + if(fShape=='a') return 3; + if(fShape=='d') return 3; // WAS 4 + if(fShape=='c') return 1; + std::cout << "Module::numberOfHgcrocs() fType = " << fType + << ", fShape = " << fShape << std::endl; + assert(false); + } + + if(fType=='M') { + if(fShape=='F') return 3; + //if(fShape=='b') return 3; + std::cout << "Module::numberOfHgcrocs() fType = " << fType + << ", fShape = " << fShape << std::endl; + assert(false); + } + + if(fType=='O') { + if(fShape=='F') return 3; + if(fShape=='b') return 3; + if(fShape=='g') return 3; + if(fShape=='a') return 2; + if(fShape=='d') return 2; + if(fShape=='c') return 1; + if(fShape=='h') return 1; + if(fShape=='i') return 1; + if(fShape=='j') return 2; + if(fShape=='y') return 1; + if(fShape=='z') return 2; + std::cout << "Module::numberOfHgcrocs() fType = " << fType + << ", fShape = " << fShape << std::endl; + assert(false); + } + + std::cout << "Module::numberOfHgcrocs() fShape = " << fShape << std::endl; + assert(false); + } + + void print() { + std::cout << "Module::print() type " << fType << ", shape " << fShape + << ", type number " << type() << ", number of HGCROCs " + << numberOfHgcrocs() + << ", number of DAQ Elinks " + << numberOfDaqElinks() + << ", number of TPG Elinks " + << numberOfTpgElinks() << std::endl; + } + + char fShape; + char fType; + int fU; + int fV; + std::vector vHgcroc; +}; + +class Motherboard { +public: + Motherboard() { + fInner=0; + fMiddle=0; + fOuter=0; + fMean=0.0; + fRms=0.0; + } + + void read(std::istream &fin) { + //fin >> fNumber >> fInner >> fMiddle >> fOuter >> fMean >> fRms; + fin >> fNumber + >> fInner >> fInnerShapes + >> fMiddle >> fMiddleShapes + >> fOuter >> fOuterShapes >> fMean >> fRms; + + print(); + assert(fin); + } + + unsigned numberOfDaqElinks() const { + unsigned l(0); + for(unsigned m(0);m36*fEconD) { + std::cout << "Layer " << fLayer << ", links = " << l << " fEconD = " << fEconD << " "; + for(unsigned m(0);m36*fEconT) { + std::cout << "Layer " << fLayer << ", links = " << l << " fEconT = " << fEconT << " "; + for(unsigned m(0);m vModule; + + double fMean; + double fRms; +}; + + +class Slink { +public: + Slink() { + } + + unsigned numberOfHgcrocs() const { + unsigned n(0); + for(unsigned m(0);m vMotherboard; +}; + + +int main(int argc, char* argv[]) { + + unsigned who(1); + std::string sWho("Langford"); + sWho+=(who==0?"0":"1"); + TFileHandler tfh(sWho); + + char x; + unsigned ly,nm,md,ninner,nmiddle,nouter; + std::string sShape; + + if(who==0) { + std::ifstream fin("HGCalOccupancy_SiOnly_V3.txt"); + + std::vector vMB[50]; + + for(unsigned l(0);l<36;l++) { + fin >> x >> ly >> nm; + assert(x=='L'); + assert(ly==l); + + for(unsigned m(0);m> x; + assert(x=='M'); + //assert(md==m+1); + + vMB[l].push_back(Motherboard()); + vMB[l][vMB[l].size()-1].fLayer=l; + vMB[l][vMB[l].size()-1].read(fin); + assert(vMB[l][vMB[l].size()-1].fNumber==m+1); + + vMB[l][vMB[l].size()-1].print(); + + + unsigned mTot(vMB[l][vMB[l].size()-1].fInner+ + vMB[l][vMB[l].size()-1].fMiddle+ + vMB[l][vMB[l].size()-1].fOuter); + for(unsigned i(0);i> x; + assert(x=='W'); + + fin >> x >> md; + assert(x=='M'); + assert(md==m+1); + + std::vector &v(vMB[l][vMB[l].size()-1].vModule); + v.push_back(Module()); + v[v.size()-1].read(fin); + } + } + } + + return 0; + // IMPROVE + for(unsigned l(36);l<50;l++) { + for(unsigned m(0);m<(l<44?7:5);m++) { + vMB[l].push_back(vMB[l-1][m]); + vMB[l][m].fLayer=l; + vMB[l][m].fMean=0.99*vMB[l-1][m].fMean; + } + } + + // Plot basic properties + TH1F *hNMB=new TH1F("NMB",";Layer;Number of motherboards",50,0,50); + TH2F *hMean=new TH2F("Mean",";Layer;Motherboard;Mean",50,0,50,25,0,25); + TH2F *hRms=new TH2F("Rms",";Layer;Motherboard;RMS",50,0,50,25,0,25); + TH2F *hHgcroc=new TH2F("Hgcroc",";Layer;Motherboard;Hgcrocs",50,0,50,25,0,25); + TH2F *hOcc=new TH2F("Occ",";Layer;Motherboard;Fractional occupancy",50,0,50,25,0,25); + TH2F *hWords=new TH2F("Words",";Layer;Motherboard;Mean words",50,0,50,25,0,25); + + unsigned nTotalMB(0); + double totalMean(0.0); + + for(unsigned l(0);l<50;l++) { + nTotalMB+=vMB[l].size(); + + hNMB->Fill(l,vMB[l].size()); + for(unsigned m(0);mFill(l,m,vMB[l][m].fMean); + hRms->Fill(l,m,vMB[l][m].fRms); + hHgcroc->Fill(l,m,vMB[l][m].numberOfHgcrocs()); + hOcc->Fill(l,m,vMB[l][m].fMean/vMB[l][m].numberOfChannels()); + hWords->Fill(l,m,vMB[l][m].meanWords()); + } + } + + std::cout << "Sector, HGCAL total number of motherboards = " << nTotalMB + << ", " << 12*nTotalMB + << ", sector, HGCAL total mean = " << totalMean + << ", " << 12.0*totalMean << std::endl; + + TH2F *hOcc2=new TH2F("Occ2",";Layer;Motherboard",50,0,50,25,0,25); + + for(unsigned l(0);l<50;l++) { + std::vector vDone(vMB[l].size()); + for(unsigned n(0);nFill(l,m,mMax/vMB[l][nMax].numberOfChannels()); + } + } + + std::vector mbSort; + + std::vector vDone2[50]; + for(unsigned l(0);l<50;l++) { + for(unsigned m(0);m0) assert(mbSort[m].meanWords()<=mbSort[m-1].meanWords()); + hSortMean->Fill(m,mbSort[m].fMean); + hSortRms->Fill(m,mbSort[m].fRms); + hSortRM->Fill(m,mbSort[m].fRms/mbSort[m].fMean); + hSortWords->Fill(m,mbSort[m].meanWords()); + hSortWordsVsMean->Fill(mbSort[m].fMean,mbSort[m].meanWords()); + } + + Slink mSlink[144]; + for(unsigned s(0);s<144;s++) { + mSlink[s].vMotherboard.push_back(mbSort[s]); + } + + for(unsigned m(144);mmSlink[s].meanWords()) { + mMin=mSlink[s].meanWords(); + sMin=s; + } + } + } + mSlink[sMin].vMotherboard.push_back(mbSort[m]); + } + + TH1F *hSlinkNum=new TH1F("SlinkNum",";Slink;Number of motherboards",144,0,144); + TH1F *hSlinkMean=new TH1F("SlinkMean",";Slink;Mean of motherboards",144,0,144); + TH1F *hSlinkRms=new TH1F("SlinkRms",";Slink;Rms of motherboards",144,0,144); + TH1F *hSlinkWords=new TH1F("SlinkWords",";Slink;Mean words of motherboards",144,0,144); + + unsigned nCheck(0); + for(unsigned s(0);s<144;s++) { + hSlinkMean->Fill(s,mSlink[s].mean()); + hSlinkRms->Fill(s,mSlink[s].rms()); + hSlinkWords->Fill(s,mSlink[s].meanWords()); + hSlinkNum->Fill(s,mSlink[s].vMotherboard.size()); + nCheck+=mSlink[s].vMotherboard.size(); + } + assert(nCheck==nTotalMB); + + //////////////////////////////////////////////////////////////// + + } else if(who==1) { + std::ifstream fin("Bloch_concat_tweaked.txt"); + assert(fin); + + std::vector vMotherboard[50]; + + char array[1024]; + unsigned a,b,c,d,e; + double f,g; + std::string s; + + unsigned nTotalMotherboards(0); + + unsigned nMinDaqLinks(0); + unsigned nTotalDaqLinks(0); + unsigned nRateEconDaqLinks(0); + unsigned nBigRateDaqLinks(0); + unsigned nBigEconDaqLinks(0); + unsigned nNoEconDaqLinks0(0); + unsigned nNoEconDaqLinks1(0); + unsigned nNoEconDaqLinks2(0); + + unsigned nTotalTpgLinks(0); + unsigned nTotalTpgLinksMin(0); + unsigned nTotalTpgLinksSTC0(0); + unsigned nTotalTpgLinksSTC1(0); + unsigned nTotalTpgLinksSTC2(0); + unsigned nTotalTpgLinksDirect0(0); + unsigned nTotalTpgLinksDirect1(0); + unsigned nTotalTpgLinksDirect2(0); + unsigned nTotalTpgLinksDirect000(0); + unsigned nTotalTpgLinksDirect100(0); + unsigned nTotalTpgLinksDirect110(0); + unsigned nTotalTpgLinksDirect111(0); + + unsigned nLayer,nLayerOld(0); + fin >> nLayer; + assert(nLayer==0); + fin >> a >> b >> c >> d; + + fin >> nLayer; + assert(nLayer==0); + + TH1F *hTotalMbVsLayer=new TH1F("TotalMbVsLayer",";Layer;Number of motherboards",36,0,36); + TH1F *hDoubleMbVsLayer=new TH1F("DoubleMbVsLayer",";Layer;Number of motherboards",36,0,36); + TH1F *hLinksVsLayer=new TH1F("LinksVsLayer",";Layer;Number of lpGBTs",36,0,36); + TH1F *hRate=new TH1F("Rate",";Rate (Gbit/s);Number of lpGBTs",80,0,8); + TH2F *hRateVsLayer=new TH2F("RateVsLayer",";Layer;Rate (Gbit/s);Number of lpGBTs",36,0,36,80,0,8); + + unsigned nLayerMBs(0); + while(fin) { + assert(nLayer<36); + fin >> a >> b >> c >> f >> d >> e >> g >> s; + assert(a==1+vMotherboard[nLayer].size()); + assert(b>0 && b<3); + assert(b==c); + assert(d>0 && d<3); + assert(e>=d && e<=3*d); + + //std::cout << s << std::endl; + + nTotalDaqLinks+=b; + if(nLayer>=28 || (nLayer%2)==0) nTotalTpgLinks+=e; + + nTotalMotherboards++; + + Motherboard motherboard; + motherboard.fLayer=nLayer; + motherboard.fEconD=b; + motherboard.fEconT=d; + motherboard.fEconTLpGbts=e; + + hTotalMbVsLayer->Fill(nLayer); + if(b==2) hDoubleMbVsLayer->Fill(nLayer); + hLinksVsLayer->Fill(nLayer,b); + if(b==1) { + hRate->Fill(f); + hRateVsLayer->Fill(nLayer,f); + } else { + hRate->Fill(f/b); + hRate->Fill(f/b); + hRateVsLayer->Fill(nLayer,f/b); + hRateVsLayer->Fill(nLayer,f/b); + } + + + unsigned nSixths0(0); + unsigned nSixths1(0); + unsigned nSixths2(0); + + unsigned nDirect000(0); + unsigned nDirect100(0); + unsigned nDirect110(0); + unsigned nDirect111(0); + unsigned nDirect2(0); + + for(unsigned i(0);i6.3 || motherboard.numberOfHgcrocs()>18?2:1); + nBigRateDaqLinks+=(f>6.3?2:1); + nBigEconDaqLinks+=(motherboard.numberOfHgcrocs()+17)/18; + nNoEconDaqLinks0+=(motherboard.numberOfHgcrocs()+13)/14; + nNoEconDaqLinks1+=(motherboard.numberOfHgcrocs()+6)/7; + nNoEconDaqLinks2+=(2*motherboard.numberOfHgcrocs()+6)/7; + + if(nLayer>=28 || (nLayer%2)==0) { + unsigned minLpGbts((motherboard.numberOfTpgElinks()+35)/36); + unsigned direct0((motherboard.numberOfHgcrocs()+13)/14); + unsigned direct1((motherboard.numberOfHgcrocs()+6)/7); + std::cout << "Number of modules = " << s.size()/2 + << ", threshold links = " << e + << ", minimum links = " << minLpGbts + << ", STC0 links = " << (nSixths0+5)/6 + << ", STC1 links = " << (nSixths1+5)/6 + << ", STC2 links = " << (nSixths2+5)/6 + << ", HGCROCs = " << motherboard.numberOfHgcrocs() + << ", direct0 links = " << direct0 + << ", direct1 links = " << direct1 + << ", direct2 links = " << (nDirect2+6)/7 << std::endl; + nTotalTpgLinksMin+=minLpGbts; + nTotalTpgLinksSTC0+=std::max((nSixths0+5)/6,minLpGbts); + nTotalTpgLinksSTC1+=std::max((nSixths1+5)/6,minLpGbts); + nTotalTpgLinksSTC2+=std::max((nSixths2+5)/6,minLpGbts); + nTotalTpgLinksDirect0+=direct0; + nTotalTpgLinksDirect1+=direct1; + nTotalTpgLinksDirect2+=(nDirect2+6)/7; + + nTotalTpgLinksDirect000+=(nDirect000+13)/14; + nTotalTpgLinksDirect100+=(nDirect100+13)/14; + nTotalTpgLinksDirect110+=(nDirect110+13)/14; + nTotalTpgLinksDirect111+=(nDirect111+13)/14; + } + + motherboard.numberOfDaqElinks(); + motherboard.numberOfTpgElinks(); + vMotherboard[nLayer].push_back(motherboard); + + fin.getline(array,1024); + //std::cout << array << std::endl; + + fin >> nLayer; + if(nLayer==nLayerOld+1) { + nMinDaqLinks+=(nLayerMBs+17)/18; + nLayerMBs=0; + + std::cout << std::endl << "Changing from layer " << nLayerOld << " to " << nLayer << std::endl << std::endl; + fin >> a >> b >> c >> d; + nLayerOld++; + fin >> nLayer; + } + } + nMinDaqLinks+=(nLayerMBs+17)/18; + + std::cout << "Total number of motherboards per sector = " << nTotalMotherboards + << ", for whole HGCAL = " << 12*nTotalMotherboards << std::endl; + + std::cout << "Min number of DAQ links per sector = " << nMinDaqLinks + << ", for whole HGCAL = " << 12*nMinDaqLinks << std::endl; + std::cout << "Total number of DAQ links per sector = " << nTotalDaqLinks + << ", for whole HGCAL = " << 12*nTotalDaqLinks << std::endl; + std::cout << "Total number of rate+ECON DAQ links per sector = " << nRateEconDaqLinks + << ", for whole HGCAL = " << 12*nRateEconDaqLinks << std::endl; + std::cout << "Total number of rate ECON DAQ links per sector = " << nBigRateDaqLinks + << ", for whole HGCAL = " << 12*nBigRateDaqLinks << std::endl; + std::cout << "Total number of big ECON DAQ links per sector = " << nBigEconDaqLinks + << ", for whole HGCAL = " << 12*nBigEconDaqLinks << std::endl; + std::cout << "Total number of no ECON0 DAQ links per sector = " << nNoEconDaqLinks0 + << ", for whole HGCAL = " << 12*nNoEconDaqLinks0 << std::endl; + std::cout << "Total number of no ECON1 DAQ links per sector = " << nNoEconDaqLinks1 + << ", for whole HGCAL = " << 12*nNoEconDaqLinks1 << std::endl; + std::cout << "Total number of no ECON2 DAQ links per sector = " << nNoEconDaqLinks2 + << ", for whole HGCAL = " << 12*nNoEconDaqLinks2 << std::endl; + std::cout << std::endl; + + std::cout << "Total number of threshold TPG links per sector = " << nTotalTpgLinks + << ", for whole HGCAL = " << 12*nTotalTpgLinks << std::endl; + std::cout << "Total number of min TPG links per sector = " << nTotalTpgLinksMin + << ", for whole HGCAL = " << 12*nTotalTpgLinksMin << std::endl; + std::cout << "Total number of STC0 TPG links per sector = " << nTotalTpgLinksSTC0 + << ", for whole HGCAL = " << 12*nTotalTpgLinksSTC0 << std::endl; + std::cout << "Total number of STC1 TPG links per sector = " + << nTotalTpgLinksSTC1 + << ", for whole HGCAL = " << 12*nTotalTpgLinksSTC1 << std::endl; + std::cout << "Total number of STC2 TPG links per sector = " + << nTotalTpgLinksSTC2 + << ", for whole HGCAL = " << 12*nTotalTpgLinksSTC2 << std::endl; + std::cout << "Total number of direct0 TPG links per sector = " + << nTotalTpgLinksDirect0 + << ", for whole HGCAL = " << 12*nTotalTpgLinksDirect0 << std::endl; + std::cout << "Total number of direct1 TPG links per sector = " + << nTotalTpgLinksDirect1 + << ", for whole HGCAL = " << 12*nTotalTpgLinksDirect1 << std::endl; + std::cout << "Total number of direct2 TPG links per sector = " + << nTotalTpgLinksDirect2 + << ", for whole HGCAL = " << 12*nTotalTpgLinksDirect2 << std::endl; + std::cout << "Total number of direct000 TPG links per sector = " + << nTotalTpgLinksDirect000 + << ", for whole HGCAL = " << 12*nTotalTpgLinksDirect000 << std::endl; + std::cout << "Total number of direct100 TPG links per sector = " + << nTotalTpgLinksDirect100 + << ", for whole HGCAL = " << 12*nTotalTpgLinksDirect100 << std::endl; + std::cout << "Total number of direct110 TPG links per sector = " + << nTotalTpgLinksDirect110 + << ", for whole HGCAL = " << 12*nTotalTpgLinksDirect110 << std::endl; + std::cout << "Total number of direct111 TPG links per sector = " + << nTotalTpgLinksDirect111 + << ", for whole HGCAL = " << 12*nTotalTpgLinksDirect111 << std::endl; + } else { + assert(false); + } + + return 0; +} diff --git a/src/SuperTCs.cpp b/src/SuperTCs.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d346414a20e5c826e682a79f9a349ad217a94be6 --- /dev/null +++ b/src/SuperTCs.cpp @@ -0,0 +1,2312 @@ +/* + g++ -Wall -std=c++0x -I. code/CalibrationTest.cc -o CalibrationTest.exe `root-config --libs --cflags` +*/ + +/* + 0 = SCs + 1 = SCs with perfect ADC + + 0 = All TCs + 0 = True energy TCs + 1 = All TCs with cell centring + 2 = Threshold TCs + 3 = Equal share (0 bits) + 4 = Equal share if > 0 (4 bits) + 5 = SuperTCs (2 bits) + + 10-17 = Fractions 1-8 bits/TC (5-26 bits) + 18 = HGCROC-like implementation of fraction 1 bit + 19 = HGCROC-like implementation of fraction 2 bits + + 20-25 = Compression 8-3 bits/TC, drop LSBs (16-32 bits) + 30-35 = Compression 8-3 bits/TC, saturate (16-32 bits) + + 40 = Fractions 1 bit/TC + 4E/5M + 41 = Fractions 2 bits/TC + 3E/3M + 42 = Fractions 2 bits/TC + 4E/2M + 43 = Fractions 2 bits/TC + 5E/1M + + + Work to limits when compressing: + SC 25 x 2^12 < 2^5 x 2^12 = 2^17 + NTC lo = 4 x 25 x 2^12 = 25 x 2^14 < 2^5 x 2^15 = 2^19 + NTC hi = 9 x 25 x 2^12 = 225 x 2^12 < 2^8 x 2^12 = 2^20 + CTC = 4 x 225 x 2^12 = 225 x 2^14 < 2^8 x 2^14 = 2^22 + DONE ! Do two TOT ranges properly + DONE ! Fix Gaussian to integrate properly + Make PU from other Gaussians plus model as density not per sensor cell + Set up various positions lo/hi radius for both densities + DONE! Change numbers of cells to work for both 2.2 and 3.3 cm sizes + Make sure we catch equal edges of fractions i.e. 1/4, 1/2, etc. correctly. + Connect up fSinTheta with fEta + +*/ + + +#include "Hack.hh" + +#include + +#include +#include +#include +#include + +#include "TFile.h" +#include "TF1.h" +#include "TH1D.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TProfile.h" +//#include "TMath.h" +#include "TE1F.hh" +#include "TH1A.hh" +#include "Random.hh" +#include "PackerFloat22.hh" +#include "TFileHandler.hh" + +const unsigned kNumberOfRadii(24); + +double cutRadius(unsigned r) { + assert(r(fHiDensity?0.5:0.05)) fPileupArray[ix][iy]=0.0; + else fPileupArray[ix][iy]=0.5+Random::random().Exp(5.0); + + fEnergyArray[ix][iy]=0.0; + } + } + } + + void increment(unsigned ix, unsigned iy, double e) { + assert(ix=4) bins=(1<<(m-3))-1; + + for(unsigned jx(2*ix);jx<2*(ix+1);jx++) { + for(unsigned jy(2*iy);jy<2*(iy+1);jy++) { + if(!(jx==mx && jy==my)) { + if(fNtcAdc[kNtcReal][jx][jy]==0) fraction=0; + else fraction=(2*bins*fNtcAdc[kNtcReal][jx][jy]-1)/sum; + + assert(fraction=4 && fraction==0) fNtcAdc[a][jx][jy]=0; + + hDeltaE1[kNtcFractions+m]->Fill(sum,(double(fNtcAdc[kNtcFractions+m][jx][jy])-double(fNtcAdc[kNtcReal][jx][jy]))/sum); + hDeltaE2[kNtcFractions+m]->Fill(sum,(double(fNtcAdc[kNtcFractions+m][jx][jy])-double(fNtcAdc[kNtcReal][jx][jy]))/sum); + } + } + } + + } else if(m==8) { + unsigned sum2=sum; + + //std::cout << "Calling 4E5M with data = " << (sum+1)/2 << std::endl; + fPacker4E5M.pack((sum+1)/2); + sum=2*fPacker4E5M.unpack(); + + // Fractions 1 bit/TC + 4E/5M + unsigned bins(2); + for(unsigned jx(2*ix);jx<2*(ix+1);jx++) { + for(unsigned jy(2*iy);jy<2*(iy+1);jy++) { + if(!(jx==mx && jy==my)) { + if(fNtcAdc[kNtcReal][jx][jy]==0) fraction=0; + else fraction=(2*bins*fNtcAdc[kNtcReal][jx][jy]-1)/sum2; + assert(fraction100000) { + std::cout << "a = " << a << ", E1 = " << fNtcAdc[kNtcReal][mx][my] << " ~ " << fNtcAdc[a][mx][my] << " = " << fNtcAdc[a][2*ix ][2*iy ] + << ", E2 = " << fNtcAdc[kNtcReal][mx][ny] << " ~ " << fNtcAdc[a][mx][ny] << " = " << fNtcAdc[a][2*ix+1][2*iy ] + << ", E3 = " << fNtcAdc[kNtcReal][nx][my] << " ~ " << fNtcAdc[a][nx][my] << " = " << fNtcAdc[a][2*ix ][2*iy+1] + << ", E4 = " << fNtcAdc[kNtcReal][nx][ny] << " ~ " << fNtcAdc[a][nx][ny] << " = " << fNtcAdc[a][2*ix+1][2*iy+1] + << ", sum = " << fNtcAdc[kNtcReal][2*ix ][2*iy ]+fNtcAdc[kNtcReal][2*ix+1][2*iy ]+fNtcAdc[kNtcReal][2*ix ][2*iy+1]+fNtcAdc[kNtcReal][2*ix+1][2*iy+1] + << " ~ " << fNtcAdc[a][2*ix ][2*iy ]+fNtcAdc[a][2*ix+1][2*iy ]+fNtcAdc[a][2*ix ][2*iy+1]+fNtcAdc[a][2*ix+1][2*iy+1] + << " ~ " << sum + << std::endl; + assert(false); + } + + if(m<8) assert((fNtcAdc[kNtcReal][2*ix][2*iy]+fNtcAdc[kNtcReal][2*ix+1][2*iy]+fNtcAdc[kNtcReal][2*ix][2*iy+1]+fNtcAdc[kNtcReal][2*ix+1][2*iy+1])==sum); + assert((fNtcAdc[a][2*ix][2*iy]+fNtcAdc[a][2*ix+1][2*iy]+fNtcAdc[a][2*ix][2*iy+1]+fNtcAdc[a][2*ix+1][2*iy+1])==sum); + } + + + void digitise(unsigned mxx) { + unsigned sigma3(unsigned(3.0*fNoise*fMip+0.5)); + //std::cout << "sigma3 = " << sigma3 << std::endl; + unsigned halfMip(unsigned(0.5*fMip+0.5)); + //std::cout << "halfMip = " << halfMip << std::endl; + unsigned nGap(300); + + // Sensor cells + + for(unsigned ix(0);ix1000) fScAdc[kScReal][ix][iy]=1000+nGap; + + } else { + + // Do two-range 9-bit TOT + fScAdc[kScReal][ix][iy]=(fScAdc[kScReal][ix][iy]+13)/25; + if(fScAdc[kScReal][ix][iy]<0x200) { + fScAdc[kScReal][ix][iy]*=25; + } else { + fScAdc[kScReal][ix][iy]=(fScAdc[kScReal][ix][iy]+2)/4; + if(fScAdc[kScReal][ix][iy]>0x1ff) fScAdc[kScReal][ix][iy]=0x1ff; + fScAdc[kScReal][ix][iy]*=100; + + //std::cout << "ix,iy = " << ix << ", " << iy << " ADC = " << fScAdc[kScReal][ix][iy] << std::endl; + } + } + + hScAdcConversion->Fill(fScAdc[kScPerfect][ix][iy],fScAdc[kScReal][ix][iy]); + hScAdcDeltaConversion->Fill(fScAdc[kScPerfect][ix][iy],int(fScAdc[kScReal][ix][iy])-int(fScAdc[kScPerfect][ix][iy])); + + // Thresholds + for(unsigned m(0);m0) nNtcTh++; + + sum+=fNtcAdc[kNtcReal][jx][jy]; + if(eMax<=fNtcAdc[kNtcReal][jx][jy]) { + eMax=fNtcAdc[kNtcReal][jx][jy]; + mx=jx; + my=jy; + } + } + } + + unsigned nx((mx%2)==0?mx+1:mx-1); + unsigned ny((my%2)==0?my+1:my-1); + + //std::cout << "Doing normal trigger cells " << mx << ", " << my << std::endl; + + fPacker4E3M.pack((sum+4)/8); + fPacker4E5M.pack((sum+1)/2); + + unsigned iNtcTh(0); + for(unsigned jx(ix);jx0) { + fNtcAdc[kNtcThresholdEqualShare][jx][jy]=(2*fPacker4E5M.unpack()+iNtcTh)/nNtcTh; + iNtcTh++; + } else fNtcAdc[kNtcThresholdEqualShare][jx][jy]=0; + } + } + assert(iNtcTh==nNtcTh); + + } + } + + //std::cout << "Done coarse trigger cells" << std::endl; + + for(unsigned ix(0);ixFill(fNtcAdc[kNtcReal][ix][iy],fNtcAdc[m][ix][iy]); + hNtcAdcConversion1[m]->Fill(fNtcAdc[kNtcReal][ix][iy],fNtcAdc[m][ix][iy]); + hNtcAdcConversion2[m]->Fill(fNtcAdc[kNtcReal][ix][iy],fNtcAdc[m][ix][iy]); + hNtcAdcDeltaConversion0[m]->Fill(fNtcAdc[kNtcReal][ix][iy],int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy])); + hNtcAdcDeltaConversion1[m]->Fill(fNtcAdc[kNtcReal][ix][iy],int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy])); + hNtcAdcDeltaConversion2[m]->Fill(fNtcAdc[kNtcReal][ix][iy],int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy])); + if(fNtcAdc[kNtcReal][ix][iy]>0) { + hNtcAdcFracConversion0[m]->Fill(fNtcAdc[kNtcReal][ix][iy],(int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy]))/(1.0*fNtcAdc[kNtcReal][ix][iy])); + hNtcAdcFracConversion1[m]->Fill(fNtcAdc[kNtcReal][ix][iy],(int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy]))/(1.0*fNtcAdc[kNtcReal][ix][iy])); + hNtcAdcFracConversion2[m]->Fill(fNtcAdc[kNtcReal][ix][iy],(int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy]))/(1.0*fNtcAdc[kNtcReal][ix][iy])); + } + } + } + //std::cout << "Ntc method " << m << " sum = " << nTcSum[m] << ", Sc sum = " << nScSum << std::endl; + if(m!=kNtcEqualShare && m!=kNtcThreshold && m!=kNtcThresholdEqualShare && m!=kNtcSuperTcs && m0) chiSq[1]+=pow((int(fNtcAdc[m][jx][jy])-int(fNtcAdc[kNtcReal][jx][jy]))/(1.0*fNtcAdc[kNtcReal][jx][jy]),2); + } + } + hNtcAdcChiSq0[m]->Fill(chiSq[0]); + hNtcAdcChiSq1[m]->Fill(chiSq[1]); + + numChiSq0[m]++; + sumChiSq0[m]+=chiSq[0]; + numChiSq1[m]++; + sumChiSq1[m]+=chiSq[1]; + } + } + + hNtcAdcChiSqAvg0->SetBinContent(m+1,sumChiSq0[m]/numChiSq0[m]); + hNtcAdcChiSqAvg1->SetBinContent(m+1,sumChiSq1[m]/numChiSq1[m]); + } + + +#ifdef CRAP + + + + unsigned nSum[5]={0,0,0,0,0}; + unsigned nSum2[kNumberOfNtcMethods]={0,0,0,0,0,0,0}; + unsigned nSum3[kNumberOfCtcMethods]={0,0,0,0}; + + + + for(unsigned ix(0);ix=mipT2) fNtcAdc[1][ix][iy]=fNtcAdc[0][ix][iy]; + } + } + + for(unsigned ix(0);ixFill(ix,iy,a.scTotalEnergy(ix,iy)); + hTrueEnergyTcCm->Fill(a.sensorCellPosition(ix),a.sensorCellPosition(iy),a.scTrueEnergy(ix,iy)/tcArea); + + h2[0][1]->Fill(ix,iy,a.scAdcEnergy(ix,iy,0)); + } + } + + for(unsigned ix(0);ixFill(ix,iy,a.tcAdcEnergy(ix,iy,m)/tcArea); + hTcAdcEnergyCm[m]->Fill(a.triggerCellPosition(ix),a.triggerCellPosition(iy),a.tcAdcEnergy(ix,iy,m)/tcArea); + } + } + } + + tcArea*=4.0; + /* + for(unsigned ix(0);ixFill(ix,iy,a.ctcAdcEnergy(ix,iy,m)/tcArea); + } + } + } + + tcArea*=(a.fHiDensity?2.0:4.0); + + for(unsigned ix(0);ixFill(ix,iy,a.htcAdcEnergy(ix,iy)/tcArea); + } + } + } + */ +} + +int main(int argc, const char **argv) { + unsigned place(0); + bool ez(true); + bool pu(true); + double fLayerZ(330.0); + + double fPt=40.0; + //fPt=1000.0; + double fPtSpread=0.5; + //double fEta=3.0; + double fMipsPerGeV(10.0); + + if(argc>4) { + std::cerr << "Usage: " << argv[0] << " h or l" << std::endl; + return 1; + } + + if(argc==2) { + std::istringstream sin(argv[1]); + sin >> place; + } + + if(argc==3) { + if(argv[2][0]=='e') ez=true; + if(argv[2][0]=='z') ez=false; + } + + if(argc==4) { + if(argv[3][0]=='p') pu=true; + if(argv[3][0]=='u') pu=false; + } + + TFileHandler tfh(std::string("SuperTCs")+(place==0?"0":"1")); + + PackerFloat22 aPacker(4,7); + + unsigned bins(1<<22); + TH1F *hPackerA=new TH1F("ScPackerA",";Input;Output",bins,0,bins); + TH1F *hPackerD=new TH1F("ScPackerD",";Input;Output-Input",bins,0,bins); + TH1F *hPackerR=new TH1F("ScPackerR",";Input;(Output-Input)/Input",bins,0,bins); + + for(unsigned i(0);iSetBinContent(i+1,op); + hPackerD->SetBinContent(i+1,int(op)-int(i)); + if(i>0) hPackerR->SetBinContent(i+1,(int(op)-int(i))/(1.0*i)); + } + + PackerFloat22 bPacker(4,3); + + bins=bPacker.dataMax()+1; + TH1F *hBPackerA=new TH1F("ScBPackerA",";Input;Output",bins,-0.5,bins-0.5); + TH1F *hBPackerD=new TH1F("ScBPackerD",";Input;Output-Input",bins,-0.5,bins-0.5); + TH1F *hBPackerR=new TH1F("ScBPackerR",";Input;(Output-Input)/Input",bins,-0.5,bins-0.5); + TH1F *hCPackerA=new TH1F("ScCPackerA",";Input;Output",bins,-0.5,bins-0.5); + TH1F *hCPackerD=new TH1F("ScCPackerD",";Input;Output-Input",bins,-0.5,bins-0.5); + TH1F *hCPackerR=new TH1F("ScCPackerR",";Input;(Output-Input)/Input",bins,-0.5,bins-0.5); + + for(unsigned i(0);iSetBinContent(i+1,op); + hBPackerD->SetBinContent(i+1,int(op)-int(i)); + if(i>0) hBPackerR->SetBinContent(i+1,(int(op)-int(i))/(1.0*i)); + + op=bPacker.unpackTimes2(); + hCPackerA->SetBinContent(i+1,op/2.0); + hCPackerD->SetBinContent(i+1,int(op)/2.0-int(i)); + if(i>0) hCPackerR->SetBinContent(i+1,(int(op)/2.0-int(i))/(1.0*i)); + } + + //const double dTriggerCells(sqrt(4.908)); + //const unsigned nSensorCells[2]={36,24}; + //const double dSensorCells[2]={dTriggerCells/3.0,dTriggerCells/2.0}; + + //const unsigned nDensity(0); + + Array theArray(place,ez,pu); + + PhotonModel thePhotonModel[3]; + thePhotonModel[0]=PhotonModel(0,0.61); + //thePhotonModel[1]=PhotonModel(1,3.0); + thePhotonModel[1]=PhotonModel(1,0.059); + thePhotonModel[2]=PhotonModel(2,1.05); + unsigned usePhotonModel(1); + + double nrg(200.0); + + TH2F *hCentre=new TH2F("Centre",";x (cm);y(cm)", + 10*theArray.numberOfTriggerCells(), + -theArray.fArrayHalfSide,theArray.fArrayHalfSide, + 10*theArray.numberOfTriggerCells(), + -theArray.fArrayHalfSide,theArray.fArrayHalfSide); + + unsigned nScEventTotal(theArray.numberOfSensorCells()*theArray.numberOfSensorCells()); + unsigned nTcEventTotal(theArray.numberOfTriggerCells()*theArray.numberOfTriggerCells()); + + TH2F *hNtcCountVsR[2]; + hNtcCountVsR[0]=new TH2F("NtcCountVsR0", + ";Cut radius #rho;Number of trigger cells", + kNumberOfRadii,0.0,0.06,300,0,300); + hNtcCountVsR[1]=new TH2F("NtcCountVsR1", + ";Cut radius #rho;Number of trigger cells", + kNumberOfRadii,0.0,0.06,300,0,300); + + for(unsigned r(0);rFill(cutRadius(r),n[0]); + + for(unsigned jx(0);jx<100;jx++) { + for(unsigned jy(0);jy<100;jy++) { + + n[1]=0; + for(unsigned ix(0);ixFill(cutRadius(r),n[1],0.0001); + } + } + } + + + TH1F *hScReadout=new TH1F("ScReadout",";Number of sensor cells",nScEventTotal,0,nScEventTotal); + + TH1F *hScEnergyTrue=new TH1F("ScEnergyTrue",";Energy (MIPs)",100,0.0,50*nrg); + TE1F *hScEnergyIntegral=new TE1F("ScEnergyIntegral",";Radius (cm);Energy (MIPs)",100,0.0,5.0,false); + TH1F *hScEnergySpectrum=new TH1F("ScEnergySpectrum",";Energy (MIPs)",120,-2.0,10.0); + + // unsigned nr(10); + TH2F *hScResponse[Array::kNumberOfScMethods][kNumberOfRadii]; + TH1A *hAScResponse[Array::kNumberOfScMethods][kNumberOfRadii]; + TH2F *hNtcResponse[2][Array::kNumberOfNtcMethods][kNumberOfRadii]; + TH1A *hANtcResponse[2][Array::kNumberOfNtcMethods][kNumberOfRadii]; + + /* + TH2F *hStcResponse[kNumberOfRadii],*hFtcResponse[kNumberOfRadii]; + TH2F *hTrueCtcResponse[kNumberOfRadii],*hCtcResponse[Array::kNumberOfCtcMethods][kNumberOfRadii]; + TH2F *hTrueHtcResponse[kNumberOfRadii],*hHtcResponse[kNumberOfRadii]; + + TH1A *hATrueScResponse[kNumberOfRadii],*hAScResponse[kNumberOfRadii]; + TH1A *hATrueTcResponse[kNumberOfRadii],*hATcResponse[Array::kNumberOfNtcMethods][kNumberOfRadii]; + TH1A *hAStcResponse[kNumberOfRadii],*hAFtcResponse[kNumberOfRadii]; + TH1A *hATrueCtcResponse[kNumberOfRadii],*hACtcResponse[Array::kNumberOfCtcMethods][kNumberOfRadii]; + TH1A *hATrueHtcResponse[kNumberOfRadii],*hAHtcResponse[kNumberOfRadii]; + */ + + for(unsigned m(0);mFill(centre[0],centre[1]); + + unsigned nScEventReadout(0); + unsigned nTcEventReadout(0); + + // Process the event + theArray.initialise(); + + double sigma(2.0); + double sigma2(2.0*sigma*sigma); + + double trueEnergy(0.0); + double thisEnergy(fMipsPerGeV*fPt*cosh(theArray.fEta)); + //double thisEnergy((0.1+0.9*Random::random().Uniform())*nrg); + + unsigned modelBins(4); + double modelBinsSize(theArray.fSensorCellSide/modelBins); + double modelBinsArea(modelBinsSize*modelBinsSize); + + for(unsigned ix(0);ixFill(dr,de); + + e+=de; + } + } + + e*=std::max(Random::random().Gaus(1.0,fPtSpread),0.0); + trueEnergy+=e; + theArray.increment(ix,iy,e); + } + } + + if(nEvt<10) std::cout << "This, true energy = " << thisEnergy << ", " << trueEnergy << " MIPs" << std::endl; + hScEnergyTrue->Fill(trueEnergy); + + theArray.digitise(0); + if(nEvt<10) std::cout << "Completed digitisation" << std::endl; + + // Analyse the event + if(nEvt<5) histEvent(nEvt,theArray); + + for(unsigned r(0);r0) nScEventReadout++; + + unsigned m=Array::kScReal; + + hScEnergySpectrum->Fill(theArray.scTotalEnergy(ix,iy)); + hScAdcSpectrum->Fill(theArray.fScAdc[m][ix][iy]); + hScAdcEnergySpectrum->Fill(theArray.scAdcEnergy(ix,iy,m)); + + hScAdcVsEnergy0->Fill(theArray.scTotalEnergy(ix,iy),theArray.scAdcEnergy(ix,iy,m)); + hScAdcVsEnergy1->Fill(theArray.scTotalEnergy(ix,iy),theArray.scAdcEnergy(ix,iy,m)); + hScAdcVsEnergy2->Fill(theArray.scTotalEnergy(ix,iy),theArray.scAdcEnergy(ix,iy,m)); + hScAdcVsEnergy3->Fill(theArray.scTotalEnergy(ix,iy),theArray.scAdcEnergy(ix,iy,m)); + + double dx(theArray.sensorCellPosition(ix)-centre[0]); + double dy(theArray.sensorCellPosition(iy)-centre[1]); + double dr(sqrt(dx*dx+dy*dy)); + + for(unsigned m(0);mFill(nScEventReadout); + + nScTotal+=nScEventTotal; + nScReadout+=nScEventReadout; + + for(unsigned m(0);mFill(trueEnergy,scResponse[m][r]); + hAScResponse[m][r]->Fill(trueEnergy,scResponse[m][r]); + } + } + + //////////////////////////////////////////////////////////////// + + for(unsigned ix(0);ix0.0) nTcEventReadout++; + /* + if(theArray.fStcAdc[ix][iy]>0) { + unsigned nx((ix%2)==0?ix+1:ix-1); + unsigned ny((iy%2)==0?iy+1:iy-1); + assert(theArray.fStcAdc[nx][iy]==0); + assert(theArray.fStcAdc[ix][ny]==0); + assert(theArray.fStcAdc[nx][ny]==0); + + hTcAdcRatio0->Fill(1.0*theArray.fTcAdc[ix][iy]/theArray.fStcAdc[ix][iy]); + hTcAdcRatio0VsStcAdc->Fill(theArray.fStcAdc[ix][iy],1.0*theArray.fTcAdc[ix][iy]/theArray.fStcAdc[ix][iy]); + + hTcAdcRatio1->Fill(1.0*theArray.fTcAdc[nx][iy]/theArray.fTcAdc[ix][iy]); + hTcAdcRatio1->Fill(1.0*theArray.fTcAdc[ix][ny]/theArray.fTcAdc[ix][iy]); + hTcAdcRatio1VsStcAdc->Fill(theArray.fStcAdc[ix][iy],1.0*theArray.fTcAdc[nx][iy]/theArray.fTcAdc[ix][iy]); + hTcAdcRatio1VsStcAdc->Fill(theArray.fStcAdc[ix][iy],1.0*theArray.fTcAdc[ix][ny]/theArray.fTcAdc[ix][iy]); + + hTcAdcRatio2->Fill(1.0*theArray.fTcAdc[nx][ny]/theArray.fTcAdc[ix][iy]); + hTcAdcRatio2VsStcAdc->Fill(theArray.fStcAdc[ix][iy],1.0*theArray.fTcAdc[nx][ny]/theArray.fTcAdc[ix][iy]); + } + */ + for(unsigned m(0);mFill(theArray.tcAdcEnergy(ix,iy,m)); + hTcAdcVsEnergy0[m]->Fill(theArray.tcTotalEnergy(ix,iy),theArray.tcAdcEnergy(ix,iy,m)); + hTcAdcVsEnergy1[m]->Fill(theArray.tcTotalEnergy(ix,iy),theArray.tcAdcEnergy(ix,iy,m)); + hTcAdcVsEnergy2[m]->Fill(theArray.tcTotalEnergy(ix,iy),theArray.tcAdcEnergy(ix,iy,m)); + hTcAdcVsEnergy3[m]->Fill(theArray.tcTotalEnergy(ix,iy),theArray.tcAdcEnergy(ix,iy,m)); + } + /* + hStcAdcEnergySpectrum->Fill(theArray.stcAdcEnergy(ix,iy)); + hStcAdcVsEnergy0->Fill(theArray.tcTotalEnergy(ix,iy),theArray.stcAdcEnergy(ix,iy)); + hStcAdcVsEnergy1->Fill(theArray.tcTotalEnergy(ix,iy),theArray.stcAdcEnergy(ix,iy)); + hStcAdcVsEnergy2->Fill(theArray.tcTotalEnergy(ix,iy),theArray.stcAdcEnergy(ix,iy)); + hStcAdcVsEnergy3->Fill(theArray.tcTotalEnergy(ix,iy),theArray.stcAdcEnergy(ix,iy)); + + hFtcAdcEnergySpectrum->Fill(theArray.ftcAdcEnergy(ix,iy)); + hFtcAdcVsEnergy0->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ftcAdcEnergy(ix,iy)); + hFtcAdcVsEnergy1->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ftcAdcEnergy(ix,iy)); + hFtcAdcVsEnergy2->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ftcAdcEnergy(ix,iy)); + hFtcAdcVsEnergy3->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ftcAdcEnergy(ix,iy)); + */ + + double dx; + double dy; + + for(unsigned ic(0);ic<2;ic++) { + if(ic==0) { + dx=theArray.triggerCellPosition(ix)-centre[0]; + dy=theArray.triggerCellPosition(iy)-centre[1]; + } else { + unsigned icx(theArray.triggerCellNumber(centre[0])); + unsigned icy(theArray.triggerCellNumber(centre[1])); + dx=theArray.triggerCellPosition(ix)-theArray.triggerCellPosition(icx); + dy=theArray.triggerCellPosition(iy)-theArray.triggerCellPosition(icy); + } + double dr(sqrt(dx*dx+dy*dy)); + + for(unsigned m(0);mFill(nTcEventReadout); + + nTcTotal+=nTcEventTotal; + nTcReadout+=nTcEventReadout; + + for(unsigned m(0);mFill(trueEnergy,ntcResponse[0][m][r]); + hANtcResponse[0][m][r]->Fill(trueEnergy,ntcResponse[0][m][r]); + hNtcResponse[1][m][r]->Fill(trueEnergy,ntcResponse[1][m][r]); + hANtcResponse[1][m][r]->Fill(trueEnergy,ntcResponse[1][m][r]); + } + } + + //////////////////////////////////////////////////////////////// + +#ifdef CRAP + for(unsigned ix(0);ixFill(theArray.ctcAdcEnergy(ix,iy,m)); + hCtcAdcVsEnergy0->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ctcAdcEnergy(ix,iy,m)); + hCtcAdcVsEnergy1->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ctcAdcEnergy(ix,iy,m)); + hCtcAdcVsEnergy2->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ctcAdcEnergy(ix,iy,m)); + hCtcAdcVsEnergy3->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ctcAdcEnergy(ix,iy,m)); + } + + double dx(theArray.coarseCellPosition(ix)-centre[0]); + double dy(theArray.coarseCellPosition(iy)-centre[1]); + double dr(sqrt(dx*dx+dy*dy)); + + for(unsigned r(0);rFill(trueEnergy,trueCtcR[r]); + hATrueCtcResponse[r]->Fill(trueEnergy,trueCtcR[r]); + + for(unsigned m(0);mFill(trueEnergy,ctcR[m][r]); + hACtcResponse[m][r]->Fill(trueEnergy,ctcR[m][r]); + } + } + + //////////////////////////////////////////////////////////////// + + for(unsigned ix(0);ixFill(trueEnergy,trueHtcR[r]); + hATrueHtcResponse[r]->Fill(trueEnergy,trueHtcR[r]); + + hHtcResponse[r]->Fill(trueEnergy,htcR[r]); + hAHtcResponse[r]->Fill(trueEnergy,htcR[r]); + } +#endif + + } // End of event loop + + //////////////////////////////////////////////////////////////// + + std::cout << "Total, read out sensor cells = " << nScTotal + << ", " << nScReadout << ", occupancy = " + << 1.0*nScReadout/nScTotal << std::endl; + std::cout << "Total, read out trigger cells = " << nTcTotal + << ", " << nTcReadout << ", occupancy = " + << 1.0*nTcReadout/nTcTotal << std::endl; + std::cout << std::endl; + + double trueRms(sqrt(xScMs[0]/nScMs[0])); + + std::cout << "Number, sumSq and RMS for true energy sensor cells = " << nScMs[0] + << ", " << xScMs[0] << ", " << sqrt(xScMs[0]/nScMs[0]) + << ", " << sqrt(xScMs[0]/nScMs[0])/trueRms << std::endl; + std::cout << "Number, sumSq and RMS for ADC energy sensor cells = " << nScMs[1] + << ", " << xScMs[1] << ", " << sqrt(xScMs[1]/nScMs[1]) + << ", " << sqrt(xScMs[1]/nScMs[1])/trueRms << std::endl; + std::cout << std::endl; + + for(unsigned m(0);mAvg(); + TH1F *hRms=hAScResponse[m][r]->Rms(); + + std::cout << "Sc method " << m << ", radius " << fLayerZ*cutRadius(r) << ", entries and integral = " + << hAvg->GetEntries() << ", " << hAvg->Integral() << std::endl; + + if(hAvg->Integral()>100) { + hAvg->Fit("pol1"); + hAvg->GetFunction("pol1")->GetParameters(parAvg); + + hScClVsR[m]->SetBinContent(r+1,parAvg[1]); + hScPuVsR[m]->SetBinContent(r+1,parAvg[0]/parAvg[1]); + + hRms->Fit(sFit.c_str()); + hRms->GetFunction(sFit.c_str())->GetParameters(parRms); + + double typical=(*(hRms->GetFunction(sFit.c_str())))(4000)/parAvg[1]; + std::cout << "Sc method " << m << ", radius " << fLayerZ*cutRadius(r) << ", typical RMS = " << typical; + + hScEsVsR[m]->SetBinContent(r+1,typical); + + hAvg->RecursiveRemove(hAvg->GetFunction("pol1")); + hRms->RecursiveRemove(hRms->GetFunction(sFit.c_str())); + + if(minRms>typical) { + minRms=typical; + minR=fLayerZ*cutRadius(r); + std::cout << " new min" << std::endl; + } else { + std::cout << std::endl; + } + } + } + std::cout << "Sc, method " << m << ", min radius " << minR << ", typical RMS = " << minRms << " <<<=======" << std::endl << std::endl; + } + + + TH1F *hNtcEsVsR[2][Array::kNumberOfNtcMethods]; + TH1F *hNtcClVsR[2][Array::kNumberOfNtcMethods]; + TH1F *hNtcPuVsR[2][Array::kNumberOfNtcMethods]; + + for(unsigned m(0);mAvg(); + TH1F *hRms=hANtcResponse[ic][m][r]->Rms(); + + std::cout << "Ntc method " << m << ", radius " << fLayerZ*cutRadius(r) << ", entries and integral = " + << hAvg->GetEntries() << ", " << hAvg->Integral() << std::endl; + + if(hAvg->Integral()>100) { + hAvg->Fit("pol1"); + hAvg->GetFunction("pol1")->GetParameters(parAvg); + + hNtcClVsR[ic][m]->SetBinContent(r+1,parAvg[1]); + hNtcPuVsR[ic][m]->SetBinContent(r+1,parAvg[0]/parAvg[1]); + + hRms->Fit(sFit.c_str()); + hRms->GetFunction(sFit.c_str())->GetParameters(parRms); + + double typical=(*(hRms->GetFunction(sFit.c_str())))(4000)/parAvg[1]; + std::cout << "Ntc method " << m << ", radius " << fLayerZ*cutRadius(r) << ", typical RMS = " << typical; + + hNtcEsVsR[ic][m]->SetBinContent(r+1,typical); + + hAvg->RecursiveRemove(hAvg->GetFunction("pol1")); + hRms->RecursiveRemove(hRms->GetFunction(sFit.c_str())); + + if(minRms>typical) { + minRms=typical; + minR=fLayerZ*cutRadius(r); + std::cout << " new min" << std::endl; + } else { + std::cout << std::endl; + } + } + } + + std::cout << "Centre method " << ic << ", Ntc method " << m << ", min radius " << minR << ", typical RMS = " << minRms << " <<<=======" << std::endl << std::endl; + } + } + + /* + for(unsigned m(0);mAvg(); + TH1F *hRms=hACtcResponse[m][r]->Rms(); + + hAvg->Fit("pol1"); + hAvg->GetFunction("pol1")->GetParameters(parAvg); + hRms->Fit(sFit.c_str()); + hRms->GetFunction(sFit.c_str())->GetParameters(parRms); + + hCtcClVsR->SetBinContent(r+1,parAvg[0]); + hCtcPuVsR->SetBinContent(r+1,parAvg[0]); + + double typical=(*(hRms->GetFunction(sFit.c_str())))(4000)/parAvg[1]; + std::cout << "Sc, radius " << fLayerZ*cutRadius(r) << ", typical RMS = " << typical; + + hCtcEsVsR->SetBinContent(r+1,typical); + hAvg->RecursiveRemove(hAvg->GetFunction("pol1")); + hRms->RecursiveRemove(hRms->GetFunction(sFit.c_str())); + + if(minRms>typical) { + minRms=typical; + minR=fLayerZ*cutRadius(r); + std::cout << " new min" << std::endl; + } else { + std::cout << std::endl; + } + } + std::cout << "Ctc method " << m << ", min radius " << minR << ", typical RMS = " << minRms << " <<<=======" << std::endl << std::endl; + } + + + minRms=1.0e12; + minR=-1.0; + + TH1F *hHtcEsVsR=new TH1F("HtcEsVsR",";Cut radius #rho;Energy resolution (MIPs)", + kNumberOfRadii,0.0075,0.0575); + TH1F *hHtcPuVsR=new TH1F("HtcPuVsR",";Cut radius #rho;Pileup (MIPs)", + kNumberOfRadii,0.0075,0.0575); + + for(unsigned r(0);rAvg(); + TH1F *hRms=hAHtcResponse[r]->Rms(); + + hAvg->Fit("pol1"); + hAvg->GetFunction("pol1")->GetParameters(parAvg); + hRms->Fit(sFit.c_str()); + hRms->GetFunction(sFit.c_str())->GetParameters(parRms); + + hHtcPuVsR->SetBinContent(r+1,parAvg[0]); + + double typical=(*(hRms->GetFunction(sFit.c_str())))(4000)/parAvg[1]; + std::cout << "Htc, radius " << fLayerZ*cutRadius(r) << ", typical RMS = " << typical; + + hHtcEsVsR->SetBinContent(r+1,typical); + + hAvg->RecursiveRemove(hAvg->GetFunction("pol1")); + hRms->RecursiveRemove(hRms->GetFunction(sFit.c_str())); + + if(minRms>typical) { + minRms=typical; + minR=fLayerZ*cutRadius(r); + std::cout << " new min" << std::endl; + } else { + std::cout << std::endl; + } + } + + std::cout << "Htc, min radius " << minR << ", typical RMS = " << minRms << " <<<=======" << std::endl << std::endl; + } + */ + } + return 0; +}