ATLAS Offline Software
TileFitter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 //************************************************************************************
6 // Filename : TileFitter.cxx
7 // Author : F. Merritt, A. Aurisano
8 // Created : March 2004
9 //
10 // DESCRIPTION
11 //
12 //*************************************************************************************
13 
16 #include <iostream>
17 #include <iomanip>
18 #include "boost/io/ios_state.hpp"
19 
20 //std::vector<double> Amp;
21 
22 // Constructors
24  : m_iConstraint(0)
25  , m_cAmp()
26  , m_pedConst(0.0)
27  , m_sigPedF2(0.0)
28  , m_sigPileF2(0.0)
29  , m_nSamples(0)
30  , m_nParameters(0)
31  , m_iConfig(0)
32  , m_vConfig()
33  , m_SPD()
34  , m_MISPD()
35  , m_errDiag()
36 {
37 }
38 
39 TileFitter::TileFitter(bool debug, int nrow, int ncol, int index, CLHEP::HepMatrix& S, int Icon)
40  : m_iConstraint(0)
41  , m_cAmp()
42  , m_pedConst(0.0)
43  , m_sigPedF2(0.0)
44  , m_sigPileF2(0.0)
45  , m_nSamples(0)
46  , m_nParameters(0)
47  , m_iConfig(0)
48  , m_vConfig()
49  , m_SPD()
50  , m_MISPD()
51  , m_errDiag()
52 {
53 
54  // Set P_const.
55  m_iConstraint = Icon;
56  m_nParameters = nrow;
57  m_nSamples = ncol;
58  m_iConfig = index;
59 
60  // if(nrow>ncol) debug=true;
61  m_pedConst = 50.0;
62  m_sigPedF2 = std::pow(1.6 / 10., 2);
63  m_sigPileF2 = std::pow(1.6 / 50., 2);
64 
65  // std::cout << "TileFitter constructor called for NP=" << NP << ", iconfig=" << iconfig << std::endl;
66  m_SPD = S;
67 
68  //Calculate the matrices needed for fitting. Save SPD and MISPD.
69  CLHEP::HepMatrix SPDT = m_SPD.T();
70  CLHEP::HepMatrix M = m_SPD * SPDT;
71 
72  // Add constraint terms to chisquare if Iconstraint>0.
73  if (m_iConstraint > 0) {
74  if (debug) {
75  std::cout << " add f*p: F2_sigPed=" << m_sigPedF2 << ", Ped_const=" << m_pedConst << ", F2_sigPile="
76  << m_sigPileF2 << std::endl;
77  }
78  M[0][0] = M[0][0] + m_sigPedF2;
79  if (m_iConstraint > 1) {
80  for (int i = 2; i < m_nParameters; i++) {
81  M[i][i] = M[i][i] + m_sigPileF2;
82  }
83  }
84  }
85 
86  int err;
87  CLHEP::HepMatrix MI = M.inverse(err);
88  m_MISPD = MI * m_SPD;
89 
90  if (m_iConstraint > 0) {
91  CLHEP::HepVector temp(m_nParameters);
92  for (int i = 0; i < 10; i++) {
93  temp[i] = MI[0][i];
94  }
95  m_cAmp = m_pedConst * m_sigPedF2 * temp;
96  if (debug) {
97  std::cout << " Camp: ";
99  }
100  }
101 
102  for (int i = 0; i < m_nParameters; i++) {
103  double err = sqrt(MI[i][i]);
104  m_errDiag.push_back(err);
105  }
106 
107  if (debug) {
108  std::cout << " SPD: ";
109  printMat(m_SPD);
110  std::cout << " SPDT: ";
111  printMat(SPDT);
112  CLHEP::HepMatrix Ident = M * MI;
113  std::cout << " M*MI: ";
114  printMat(Ident);
115  std::cout << " MISPD: ";
116  printMat(m_MISPD);
117  std::cout << " errDiag=";
118  for (int i = 0; i < m_nParameters; i++) {
119  std::cout << m_errDiag[i] << " ";
120  }
121  std::cout << std::endl;
122  }
123 }
124 // ==============================================================================
125 int TileFitter::fitAmp(TileFilterResult &tResult, bool lDebug) {
126  int iret = -1;
127  bool debugFitAmp = lDebug;
128  // Get reference to variables which will store results.
129  double sigDig = tResult.getSigDig();
130  CLHEP::HepVector& digits = tResult.getDigRef();
131  // int nrow = digits.num_row();
132  CLHEP::HepVector& Amp = tResult.getParamRef();
133  CLHEP::HepVector& Err = tResult.getErrRef();
134  CLHEP::HepVector& residuals = tResult.getResidRef();
135  double& chisq = tResult.getChi2Ref();
136  if (debugFitAmp) {
137  std::cout << " Enter FitAmp. NP=" << m_nParameters << ", ND=" << m_nSamples << std::endl;
138  std::cout << " In FitAmp: MISPD=" << std::endl;
139  printMat(m_MISPD);
140  }
141 
142  // Multiply MISPD into digits to get fitted parameters.
143 
144  Amp = m_MISPD * digits;
145  if (m_iConstraint > 0) Amp = Amp + m_cAmp;
146 
147  // Right-multiply Fit into SPD to get predicted digits.
148  CLHEP::HepVector predict(m_nSamples);
149  predict = m_SPD.T() * Amp;
150 
151  // Subtract digits from predict to get residuals, sum to get chisq*sig.
152  residuals = digits - predict;
153  chisq = residuals.normsq() / (sigDig * sigDig);
154  if (m_iConstraint > 0) chisq += std::pow(Amp[0] - m_pedConst, 2) * m_sigPedF2;
155  if (m_iConstraint > 1) {
156  for (int iamp = 2; iamp < m_nParameters; iamp++) {
157  chisq = chisq + Amp[iamp] * Amp[iamp] * m_sigPileF2;
158  }
159  }
160 
161  // Get error on each parameter from ErrDiag.
162  CLHEP::HepVector newErr(m_nParameters);
163  for (int ip = 0; ip < m_nParameters; ip++) {
164  newErr[ip] = sigDig * m_errDiag[ip];
165  }
166  Err = newErr;
167 
168  if (debugFitAmp) {
169  std::cout << " predict: ";
170  printVec(predict);
171  std::cout << " residuals: ";
172  printVec(residuals);
173  //Check chisq.
174  // double dig2 = digits.normsq();
175  // double chisq2 = dig2 - dot(M*Amp,Amp);
176  // std::cout << " chisq=" << chisq << ", chisq2=" << chisq2 << std::endl;
177  }
178  iret = 0;
179  return iret;
180 }
181 
182 // ==============================================================================
184  return m_iConfig;
185 }
186 // ==============================================================================
187 std::vector<double>& TileFitter::getErr() {
188  return m_errDiag;
189 }
190 // ==============================================================================
191 void TileFitter::printMat(CLHEP::HepMatrix &mat) {
192  boost::io::ios_base_all_saver coutsave(std::cout);
193  int nrow = mat.num_row();
194  int ncol = mat.num_col();
195  std::cout << " nrow=" << nrow << ", ncol=" << ncol << std::endl;
196  std::streamsize oldprec = std::cout.precision(4);
197  for (int irow = 0; irow < nrow; irow++) {
198  std::cout << " irow=" << irow << ": Mat=";
199  for (int icol = 0; icol < ncol; icol++) {
200  std::cout << std::setw(7) /* << std::setprecision(4) */<< mat[irow][icol];
201  }
202  std::cout << std::endl;
203  }
204  std::cout.precision(oldprec);
205 }
206 // ==============================================================================
207 void TileFitter::printVec(CLHEP::HepVector &vec) {
208  int nr = vec.num_row();
209  int nc = vec.num_col();
210  std::cout << " nr=" << nr << ", nc=" << nc << ", vec=";
211  for (int i = 0; i < nr; i++) {
212  std::cout << vec[i] << " ";
213  }
214  std::cout << std::endl;
215 }
TileFitter::m_SPD
CLHEP::HepMatrix m_SPD
Definition: TileFitter.h:60
TileFitter::m_sigPileF2
double m_sigPileF2
Definition: TileFitter.h:54
TileFitter::m_pedConst
double m_pedConst
Definition: TileFitter.h:52
TileFilterResult::getChi2Ref
double & getChi2Ref()
Definition: TileFilterResult.cxx:74
TileFitter::m_MISPD
CLHEP::HepMatrix m_MISPD
Definition: TileFitter.h:63
index
Definition: index.py:1
mat
GeoMaterial * mat
Definition: LArDetectorConstructionTBEC.cxx:55
TileFitter.h
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
TileFilterResult::getResidRef
CLHEP::HepVector & getResidRef()
Definition: TileFilterResult.cxx:70
TileFilterResult::getParamRef
CLHEP::HepVector & getParamRef()
Definition: TileFilterResult.cxx:62
TileFitter::m_errDiag
std::vector< double > m_errDiag
Definition: TileFitter.h:64
vec
std::vector< size_t > vec
Definition: CombinationsGeneratorTest.cxx:12
JetTiledMap::S
@ S
Definition: TiledEtaPhiMap.h:44
atlasStyleMacro.icol
int icol
Definition: atlasStyleMacro.py:13
TileFitter::m_iConstraint
int m_iConstraint
Definition: TileFitter.h:50
TileFilterResult::getErrRef
CLHEP::HepVector & getErrRef()
Definition: TileFilterResult.cxx:66
dqt_zlumi_pandas.err
err
Definition: dqt_zlumi_pandas.py:182
lumiFormat.i
int i
Definition: lumiFormat.py:85
TileFitter::TileFitter
TileFitter()
Definition: TileFitter.cxx:23
TileFitter::printMat
void printMat(CLHEP::HepMatrix &mat)
Definition: TileFitter.cxx:191
TileFitter::getIndex
int getIndex()
Definition: TileFitter.cxx:183
find_tgc_unfilled_channelids.ip
ip
Definition: find_tgc_unfilled_channelids.py:3
TileFitter::getErr
std::vector< double > & getErr()
Definition: TileFitter.cxx:187
TileFitter::m_iConfig
int m_iConfig
Definition: TileFitter.h:58
debug
const bool debug
Definition: MakeUncertaintyPlots.cxx:53
TileFitter::m_sigPedF2
double m_sigPedF2
Definition: TileFitter.h:53
TileFilterResult.h
TileFitter::fitAmp
int fitAmp(TileFilterResult &tResult, bool lDebug=false)
Definition: TileFitter.cxx:125
DeMoScan.index
string index
Definition: DeMoScan.py:364
TileFilterResult::getDigRef
CLHEP::HepVector & getDigRef()
Definition: TileFilterResult.cxx:46
TileFitter::m_nSamples
int m_nSamples
Definition: TileFitter.h:56
TileFilterResult
Auxiliary class for TileRawChannelMakerManyAmps.
Definition: TileFilterResult.h:22
TileFitter::m_nParameters
int m_nParameters
Definition: TileFitter.h:57
TileFitter::printVec
void printVec(CLHEP::HepVector &vec)
Definition: TileFitter.cxx:207
plotBeamSpotMon.nc
int nc
Definition: plotBeamSpotMon.py:83
TileFilterResult::getSigDig
double getSigDig() const
Definition: TileFilterResult.cxx:41
TileFitter::m_cAmp
CLHEP::HepVector m_cAmp
Definition: TileFitter.h:51