ATLAS Offline Software
TileFilterTester.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // **************************************************************************************************
6 // Filename: TileFilterTester.cxx
7 // Author: F. Merritt
8 // Created: April 2004
9 
10 // DESCRIPTION
11 // This class is to be called from the TileHitToRawChannel initialization section,
12 // after TileFilterManager has been instantiated. A reference to TileFilterManager
13 // is passed in the call, so it is possible to thoroughly test the program outside
14 // of the event loop.
15 //
16 // *************************************************************************************************
17 #include "GaudiKernel/Chrono.h"
18 //#include "GaudiKernel/ISvcLocator.h"
19 
20 #include "TileRecUtils/TileFilterTester.h"
21 #include "TileRecUtils/TileFilterManager.h"
22 #include "TileRecUtils/TileFitter.h"
23 #include "TileRecUtils/TileFilterResult.h"
24 //#include "TileConditions/TileInfo.h"
25 
26 #include <CLHEP/Random/Randomize.h>
27 #include <CLHEP/Matrix/Matrix.h>
28 #include <CLHEP/Matrix/Vector.h>
29 
30 using namespace CLHEP;
31 
32 #include <iostream>
33 #include <iomanip>
34 
35 // ==============================================================================================
36 
37 // Constructor
38 TileFilterTester::TileFilterTester(TileFilterManager* filterManager, int filterMode, int filterTest, bool debug)
39  : m_debug(debug)
40  , m_filterMode(filterMode)
41  , m_filterTest(filterTest)
42  , m_cMode(2)
43  , m_nCross(9)
44  , m_tileFilterManager(filterManager)
45 {
46 
47  std::cout << " Enter TileFilterTester constructor. April" << std::endl;
48 
49  m_nPileup = m_filterTest - 2;
50  if (m_nPileup < 0) m_nPileup = 0;
51  if (m_nPileup > m_nCross) m_nPileup = m_nCross;
52  m_nAmp = m_nPileup + 1;
53  m_iAmpVec.reserve(m_nAmp);
54  m_ampVec.reserve(m_nAmp);
55  m_nConfig = m_tileFilterManager->getNfitIndex();
56  if (m_debug) {
57  std::cout << " TileFilterTester: Cmode=" << m_cMode << ", Npileup=" << m_nPileup << std::endl;
58  for (int i = 0; i < m_nAmp; i++) {
59  std::cout << " Namp=" << i << " => Nconfig=" << m_nConfig[i] << std::endl;
60  }
61  }
62 
63  return;
64 }
65 
66 // Destructor
67 TileFilterTester::~TileFilterTester() {
68 }
69 
70 //
71 // **********************************************************************************************
72 //
73 // Generate events for testing TileFilterManager and related classes.
74 void TileFilterTester::genEvents(int nEvent) {
75  //
76  // Set up the vector of Hit amplitudes that will be used to generate fake events.
77  // This are indexed according to the "crossing index", 0=intime and 1 to Ncross-1
78  // for the pileup amplitudes. All amps are in RC units.
79  //
80  // Define the range of amplitudes to be generated.
81  double ApileMin = 30.;
82  double ApileMax = 400.;
83  double AinTime = 100.;
84  double digSigma = 1.6;
85  // Define the pileup configuration that will be used to generate events.
86  if (m_cMode == 1) { // Fixed configuration. Must define iAmpVec for first Namp elements.
87  m_iAmpVec[0] = 0;
88  m_iAmpVec[1] = 2;
89  m_iAmpVec[2] = 5;
90  }
91  bool FSM = true;
92  if (FSM) {
93  nEvent = 1000000;
94  }
95  // Initialization for sums.
96  int Ngen = 0;
97  int Nrec = 0;
98  const int ncodeSize = 8;
99  int ncode[ncodeSize];
100  int ncodeg[ncodeSize];
101  int ncodeb[ncodeSize];
102  for (int i = 0; i < ncodeSize; i++) {
103  ncode[i] = 0;
104  ncodeg[i] = 0;
105  ncodeb[i] = 0;
106  }
107  double Dsum = 0.;
108  double D2sum = 0.;
109  double E2sum = 0.;
110  int Nrecg = 0;
111  double Dsumg = 0.;
112  double D2sumg = 0.;
113  double E2sumg = 0.;
114  int Nrecb = 0;
115  double D2sumb = 0.;
116  double E2sumb = 0.;
117  int Nover = 0;
118  int Nunder = 0;
119  int Nmixed = 0;
120  int NampFinal[12];
121 
122  bool lPrint = m_debug;
123  lPrint = true;
124  bool lPrint2 = true;
125  int nPrint = 0;
126  int nPrint2 = 0;
127  for (int i = 0; i < 12; i++) {
128  NampFinal[i] = 0;
129  }
130  //
131  // ***** Start event loop. *****
132  // std::cout << " Start event loop in TileFilterTester: Nevent =" << Nevent << std::endl;
133  for (int ievent = 0; ievent < nEvent; ievent++) {
134  if (ievent > 5) m_debug = false;
135  lPrint = m_debug;
136  if (ievent == 9) lPrint = true;
137  if (ievent == 23) lPrint = true;
138  if (ievent == 235) lPrint = true;
139  if (ievent == 372) lPrint = true;
140  if (ievent == 460) lPrint = true;
141  if (ievent == 483) lPrint = true;
142  if (ievent == 501) lPrint = true;
143  if (ievent == 823) lPrint = true;
144  if (ievent == 959) lPrint = true;
145  if (ievent == 1220) lPrint = true;
146  m_iAmpVec.clear();
147  m_iAmpVec.reserve(m_nAmp);
148  int iConfig = 0;
149 
150  if (lPrint)
151  std::cout << " TileFilterTester: Start new event; ievent=" << ievent << ", Fmode=" << m_filterMode << ", Cmode="
152  << m_cMode << std::endl;
153  // Generate the pileup configuration.
154  if (m_cMode == 2) {
155  int Nindex = m_nConfig[m_nPileup];
156  double Ranflat = RandFlat::shoot();
157  iConfig = (int) (Ranflat * Nindex);
158  if (m_nPileup == 0) iConfig = 0; // Special case!
159  m_tileFilterManager->getVcross(m_nPileup, iConfig, m_iAmpVec);
160  }
161  //Generate the amplitudes.
162  m_ampVec[0] = AinTime;
163  if (m_nPileup > 0) {
164  for (int ipileup = 0; ipileup < m_nPileup; ipileup++) {
165  double Ranflat = RandFlat::shoot();
166  m_ampVec[ipileup + 1] = ApileMin + Ranflat * (ApileMax - ApileMin);
167  }
168  }
169 
170  if (lPrint) {
171  for (int iAmp = 0; iAmp < m_nAmp; iAmp++) {
172  std::cout << " i=" << iAmp << ", ipar=" << m_iAmpVec[iAmp] << ", Amp=" << m_ampVec[iAmp] << std::endl;
173  }
174  }
175  //Generate the TileDigits amplitudes.
176  int Nparam = m_nPileup + 2;
177  const int Ndig = 9;
178  HepVector digitsHep;
179  std::vector<float> digits(Ndig);
180  HepVector Param(Nparam);
181  HepMatrix SPD(Nparam, Ndig);
182  Param[0] = 50.;
183  double chisqGen = 0.;
184  for (int i = 1; i < Nparam; i++) {
185  Param[i] = m_ampVec[i - 1];
186  }
187  //int iret;
188  /*iret =*/m_tileFilterManager->makeSPD(false, m_iAmpVec, SPD);
189  HepMatrix SDP = SPD.T();
190 
191  digitsHep = SDP * Param;
192 
193  // if(lPrint) std::cout << " digits=";
194  // Convert to digits and add noise
195  double sigmaGen = 1.6;
196  for (int idig = 0; idig < Ndig; idig++) {
197  double rang = RandGauss::shoot();
198  double noise = rang * sigmaGen;
199  chisqGen += rang * rang;
200  digits[idig] = digitsHep[idig] + noise;
201  // if(lPrint) std::cout << " " << std::setw(6) << std::setprecision(2) << digitsHep[idig];
202  }
203  // if(lPrint) std::cout << std::endl;
204 
205  // *************** Now test the fitting code! ********************
206  TileFilterResult tResult(digits, digSigma);
207  int icode = m_tileFilterManager->fitDigits(tResult, lPrint);
208  if (lPrint) std::cout << "TileFilterTester.GenEvents: ievent=" << ievent << ", icode=" << icode << std::endl;
209 
210  // Compare reconstruction to generation
211  std::vector<int>& vcross = tResult.getVcrossRef();
212  int ncrgen = m_iAmpVec.size();
213  int ncrrec = vcross.size();
214  NampFinal[ncrrec] += 1;
215  int& iFitIndex = tResult.getFitIndexRef();
216  // Get the theoretical error on the in-time amplitude.
217  std::vector<double>& fitterErr = m_tileFilterManager->getFitterErr(Nparam, iConfig);
218  double Qerr = sigmaGen * fitterErr[1];
219 
220  // Define lconfigOK if reconstructed configuration = generated configuration.
221  bool lconfigOK = (ncrrec == ncrgen);
222  if (iConfig != iFitIndex) lconfigOK = false;
223  if (lconfigOK) {
224  for (int icr = 0; icr < ncrgen; icr++) {
225  if (m_iAmpVec[icr] != vcross[icr]) lconfigOK = false;
226  }
227  }
228  lPrint2 = false;
229  if (!lconfigOK) {
230  lPrint2 = true;
231  nPrint2 += 1;
232  if (nPrint2 > 100) lPrint2 = false;
233  if (lPrint2) {
234  std::cout << std::endl;
235  std::cout << " Discrepancy: ievent=" << ievent << ", ncrgen=" << ncrgen << ", ncrrec=" << ncrrec << ", icode="
236  << icode << " iConfig=" << iConfig << ", iFitIndex=" << iFitIndex << std::endl;
237  std::cout << " icrGen=";
238  for (int i = 0; i < ncrgen; i++) {
239  std::cout << " " << m_iAmpVec[i];
240  }
241  std::cout << " icrRec=";
242  for (int i = 0; i < ncrrec; i++) {
243  std::cout << " " << vcross[i];
244  }
245  double& chisq = tResult.getChi2Ref();
246  std::cout << " (chi2=" << chisq << ", chi2G=" << chisqGen << ")" << std::endl;
247 
248  for (int iAmp = 0; iAmp < m_nAmp; iAmp++) {
249  std::cout << " i=" << iAmp << ", ipar=" << m_iAmpVec[iAmp] << ", Amp=" << m_ampVec[iAmp] << std::endl;
250  }
251 
252  nPrint += 1;
253  std::cout << " digits=";
254  for (int idig = 0; idig < Ndig; idig++) {
255  std::cout << " " << std::setw(6) << std::setprecision(2) << digits[idig];
256  }
257  std::cout << std::endl;
258  tResult.printFitParam();
259  tResult.snapShot(2);
260  std::cout << std::endl;
261  if (nPrint > 10) lPrint = false;
262  }
263  }
264  // Now print out the results for comparison with generated amplitudes.
265  if (lPrint) {
266 
267  tResult.printFitParam();
268  if (nPrint2 > 20) lPrint = false;
269  }
270  double& chisq = tResult.getChi2Ref();
271  // std::vector<int>& vcross = tResult.getVcrossRef();
272  HepVector& fitParam = tResult.getParamRef();
273  HepVector& fitErr = tResult.getErrRef();
274  double diff_ch = fitParam[1] - m_ampVec[0];
275  int Npar = fitParam.num_row();
276  Ngen = Ngen + 1;
277  // double err = fitErr[1];
278  double err = Qerr;
279  Nrec += 1;
280  if (icode >= 0) ncode[icode] += 1;
281  Dsum += diff_ch;
282  D2sum += diff_ch * diff_ch;
283  E2sum += err * err;
284 
285  if (lconfigOK) {
286  Nrecg += 1;
287  if (icode >= 0) ncodeg[icode] += 1;
288  Dsumg += diff_ch;
289  D2sumg += diff_ch * diff_ch;
290  E2sumg += err * err;
291  } else {
292  Nrecb += 1;
293  if (ncrrec > ncrgen) Nover += 1;
294  if (ncrrec < ncrgen) Nunder += 1;
295  if (ncrrec == ncrgen) Nmixed += 1;
296  if (icode >= 0) ncodeb[icode] += 1;
297  D2sumb += diff_ch * diff_ch;
298  E2sumb += err * err;
299  } // end of ievent loop.
300 
301  if (lPrint)
302  std::cout << "TileFilterTester event: Npar=" << Npar << ", diff_ch =" << std::setw(6) << std::setprecision(2)
303  << diff_ch << " +-" << fitErr[1] << ", chi2=" << chisq << ", chi2Gen=" << chisqGen << std::endl << std::endl;
304  }
305  // Calculate the mean displacement and sigma of the reconstructed events.
306  std::cout << std::endl;
307  std::cout << " *** TileFilterTester Summary: Fmode=" << m_filterMode << ", Ftest=" << m_filterTest << ", NparamGen ="
308  << m_nPileup + 2 << ", Cmode=" << m_cMode << ", Nevent=" << nEvent << std::endl;
309  double rchisqCut;
310  double chiCut;
311  m_tileFilterManager->getCuts(rchisqCut, chiCut);
312  std::cout << " Cuts applied: rchisqCut=" << rchisqCut << ", chiCut=" << chiCut << std::endl;
313  std::cout << " ApileMin=" << ApileMin << ", ApileMax=" << ApileMax << ", AmpInTime=" << AinTime << std::endl;
314  std::cout << " Compare difference (rec-gen) for: Ngen=" << Ngen << ", Nrec=" << Nrec << std::endl;
315 
316  int den = Nrec ? Nrec : 1;
317  double rm = Dsum / den;
318  double errsig = pow(E2sum / den, 0.5);
319  double rsig = pow(D2sum / den, 0.5);
320 
321  std::cout << " All configurations: N=" << std::setw(7) << Nrec << ", Diff =" << std::setw(6)
322  << std::setprecision(3) << rm << ", sig =" << rsig << " (errsig=" << errsig << ")" << std::endl;
323 
324  if (Nrecg) {
325  double rmg = Dsumg / Nrecg;
326  double errsigg = pow(E2sumg / Nrecg, 0.5);
327  double rsigg = pow(D2sumg / Nrecg, 0.5);
328 
329  std::cout << " Good configurations: N=" << std::setw(7) << Nrecg << ", Diff =" << std::setw(6)
330  << std::setprecision(3) << rmg << ", sig =" << rsigg << " (errsig=" << errsigg << ")" << std::endl;
331  }
332 
333  if (Nrecb) {
334  double rmb = Dsumg / Nrecb;
335  double errsigb = pow(E2sumb / Nrecb, 0.5);
336  double rsigb = pow(D2sumb / Nrecb, 0.5);
337 
338  //if(Nrecb>3)
339  std::cout << " Bad configurations: N=" << std::setw(7) << Nrecb << ", Diff =" << std::setw(6)
340  << std::setprecision(3) << rmb << ", sig =" << rsigb << " (errsig=" << errsigb << ")" << std::endl;
341  }
342 
343  std::cout << " Nover=" << Nover << ", Nunder=" << Nunder << ", Nmixed=" << Nmixed << std::endl;
344 
345  std::cout << std::endl;
346  // Print out summary of icode counts:
347  for (int i = 0; i < ncodeSize; i++) {
348  std::cout << " icode=" << i << " ==> ncnt=" << std::setw(7) << ncode[i] << ", ncntg=" << ncodeg[i]
349  << ", ncntb=" << ncodeb[i] << std::endl;
350  }
351  std::cout << std::endl;
352  // Print out number of final amplitudes.
353  std::cout << " Number of final amplitudes:" << std::endl;
354  for (int i = 0; i < 12; i++) {
355  if (NampFinal[i] != 0) std::cout << " Namp=" << i << " =>" << NampFinal[i] << " events" << std::endl;
356  }
357 
358  return;
359 }
TileFilterManager::getCuts
void getCuts(double &rchisq, double &chi)
Definition: TileFilterManager.cxx:566
TileFilterTester.h
TileFilterTester::m_tileFilterManager
TileFilterManager * m_tileFilterManager
Definition: TileFilterTester.h:48
TileFilterTester::m_debug
bool m_debug
Definition: TileFilterTester.h:37
TileFilterResult::getChi2Ref
double & getChi2Ref()
Definition: TileFilterResult.cxx:74
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
TileFilterTester::m_nConfig
std::vector< int > m_nConfig
Definition: TileFilterTester.h:44
TileFilterManager.h
TileFilterManager::fitDigits
int fitDigits(TileFilterResult &tRes, bool lDebug=false)
Definition: TileFilterManager.cxx:160
TileFitter.h
TileFilterTester::TileFilterTester
TileFilterTester(TileFilterManager *filterManager, int filterMode, int filterTest, bool debug=false)
Definition: TileFilterTester.cxx:38
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
TileFilterResult::getParamRef
CLHEP::HepVector & getParamRef()
Definition: TileFilterResult.cxx:62
TileFilterManager
Auxiliary class for TileRawChannelMakerManyAmps.
Definition: TileFilterManager.h:34
TileFilterTester::m_nAmp
int m_nAmp
Definition: TileFilterTester.h:43
TileFilterManager::makeSPD
int makeSPD(bool debug, std::vector< int > &vcross, CLHEP::HepMatrix &SPD)
Definition: TileFilterManager.cxx:490
TileFilterResult::getErrRef
CLHEP::HepVector & getErrRef()
Definition: TileFilterResult.cxx:66
dqt_zlumi_pandas.err
err
Definition: dqt_zlumi_pandas.py:193
TileFilterTester::m_filterTest
int m_filterTest
Definition: TileFilterTester.h:39
lumiFormat.i
int i
Definition: lumiFormat.py:92
TileFilterTester::m_ampVec
std::vector< double > m_ampVec
Definition: TileFilterTester.h:46
CLHEP
STD'S.
Definition: IAtRndmGenSvc.h:19
TileFilterResult::getFitIndexRef
int & getFitIndexRef()
Definition: TileFilterResult.cxx:54
TileFilterResult::snapShot
void snapShot(int imode)
Definition: TileFilterResult.cxx:103
TileFilterTester::m_cMode
int m_cMode
Definition: TileFilterTester.h:40
TileFilterManager::getFitterErr
std::vector< double > & getFitterErr(int nParam, int iconfig)
Definition: TileFilterManager.cxx:630
TileFilterTester::m_iAmpVec
std::vector< int > m_iAmpVec
Definition: TileFilterTester.h:45
debug
const bool debug
Definition: MakeUncertaintyPlots.cxx:53
TileFilterManager::getVcross
void getVcross(int nPile, int iconfig, std::vector< int > &vcross)
Definition: TileFilterManager.cxx:574
TileFilterResult.h
TileFilterTester::~TileFilterTester
~TileFilterTester()
Definition: TileFilterTester.cxx:67
TileFilterTester::m_nPileup
int m_nPileup
Definition: TileFilterTester.h:42
TileFilterManager::getNfitIndex
std::vector< int > & getNfitIndex()
Definition: TileFilterManager.cxx:560
TileFilterTester::m_filterMode
int m_filterMode
Definition: TileFilterTester.h:38
TileFilterResult
Auxiliary class for TileRawChannelMakerManyAmps.
Definition: TileFilterResult.h:22
TileFilterTester::m_nCross
int m_nCross
Definition: TileFilterTester.h:41
TileFilterTester::genEvents
void genEvents(int nEvent)
Definition: TileFilterTester.cxx:74
WriteCellNoiseToCool.noise
noise
Definition: WriteCellNoiseToCool.py:380
TileFilterResult::getVcrossRef
std::vector< int > & getVcrossRef()
Definition: TileFilterResult.cxx:50
TileFilterResult::printFitParam
void printFitParam()
Definition: TileFilterResult.cxx:78
Generated on Sun Jun 30 2024 21:28:56 for ATLAS Offline Software by   doxygen 1.8.18