CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
SiPixelTemplate.cc
Go to the documentation of this file.
1 //
2 // SiPixelTemplate.cc Version 10.24
3 //
4 // Add goodness-of-fit info and spare entries to templates, version number in template header, more error checking
5 // Add correction for (Q_F-Q_L)/(Q_F+Q_L) bias
6 // Add cot(beta) reflection to reduce y-entries and more sophisticated x-interpolation
7 // Fix small index searching bug in interpolate method
8 // Change interpolation indexing to avoid complier complaining about possible un-initialized variables
9 // Replace containers with static arrays in calls to ysigma2 and xsigma2
10 // Add external threshold to calls to ysigma2 and xsigma2, fix parameter signal max for xsigma2
11 // Return to 5 pixel spanning but adjust boundaries to use only when needed
12 // Implement improved (faster) chi2min search that depends on pixel types
13 // Fill template arrays in single calls to this object
14 // Add qmin to the template
15 // Add qscale to match charge scales
16 // Small improvement to x-chisquare interpolation
17 // Enlarge SiPixelTemplateStore to accommodate larger templates and increased alpha acceptance (reduce PT threshold to ~200 MeV)
18 // Change output streams to conform to CMSSW info and error logging.
19 // Store x and y cluster sizes in fractional pixels to facilitate cluster splitting
20 // Add methods to return 3-d templates needed for cluster splitting
21 // Keep interpolated central 9 template bins in private storage and expand/shift in the getter functions (faster for speed=2/3) and easier to build 3d templates
22 // Store error and bias information for the simple chi^2 min position analysis (no interpolation or Q_{FB} corrections) to use in cluster splitting
23 // To save time, the gaussian centers and sigma are not interpolated right now (they aren't currently used). They can be restored by un-commenting lines in the interpolate method.
24 // Add a new method to calculate qbin for input cotbeta and cluster charge. To be used for error estimation of merged clusters in PixelCPEGeneric.
25 // Improve the charge estimation for larger cot(alpha) tracks
26 // Change interpolate method to return false boolean if track angles are outside of range
27 // Add template info and method for truncation information
28 // Change to allow template sizes to be changed at compile time
29 // Fix bug in track angle checking
30 // Accommodate Dave's new DB pushfile which overloads the old method (file input)
31 // Add CPEGeneric error information and expand qbin method to access useful info for PixelCPEGeneric
32 // Fix large cot(alpha) bug in qmin interpolation
33 // Add second qmin to allow a qbin=5 state
34 // Use interpolated chi^2 info for one-pixel clusters
35 // Fix DB pushfile version number checking bug.
36 // Remove assert from qbin method
37 // Replace asserts with exceptions in CMSSW
38 // Change calling sequence to interpolate method to handle cot(beta)<0 for FPix cosmics
39 // Add getter for pixelav Lorentz width estimates to qbin method
40 // Add check on template size to interpolate and qbin methods
41 // Add qbin population information, charge distribution information
42 //
43 //
44 // V7.00 - Decouple BPix and FPix information into separate templates
45 // Add methods to facilitate improved cluster splitting
46 // Fix small charge scaling bug (affects FPix only)
47 // Change y-slice used for the x-template to be closer to the actual cotalpha-cotbeta point
48 // (there is some weak breakdown of x-y factorization in the FPix after irradiation)
49 //
50 //
51 // V8.00 - Add method to calculate a simple 2D template
52 // Reorganize the interpolate method to extract header info only once per ID
53 // V8.01 - Improve simple template normalization
54 // V8.05 - Change qbin normalization to work better after irradiation
55 // V8.10 - Add Vavilov distribution interpolation
56 // V8.11 - Renormalize the x-templates for Guofan's cluster size calculation
57 // V8.12 - Technical fix to qavg issue.
58 // V8.13 - Fix qbin and fastsim interpolaters to avoid changing class variables
59 // V8.20 - Add methods to identify the central pixels in the x- and y-templates (to help align templates with clusters in radiation damaged detectors)
60 // Rename class variables from pxxxx (private xxxx) to xxxx_ to follow standard convention.
61 // Add compiler option to store the template entries in BOOST multiarrays of structs instead of simple c arrays
62 // (allows dynamic resizing template storage and has bounds checking but costs ~10% in execution time).
63 // V8.21 - Add new qbin method to use in cluster splitting
64 // V8.23 - Replace chi-min position errors with merged cluster chi2 probability info
65 // V8.25 - Incorporate VI's speed changes into the current version
66 // V8.26 - Modify the Vavilov lookups to work better with the FPix (offset y-templates)
67 // V8.30 - Change the splitting template generation and access to improve speed and eliminate triple index boost::multiarray
68 // V8.31 - Add correction factor: measured/true charge
69 // V8.31 - Fix version number bug in db object I/O (pushfile)
70 // V8.32 - Check for illegal qmin during loading
71 // V8.33 - Fix small type conversion warnings
72 // V8.40 - Incorporate V.I. optimizations
73 // V9.00 - Expand header to include multi and single dcol thresholds, LA biases, and (variable) Qbin definitions
74 // V9.01 - Protect against negative error squared
75 // V10.00 - Update to work with Phase 1 FPix. Fix some code problems introduced by other maintainers.
76 // V10.01 - Fix initialization style as suggested by S. Krutelyov
77 // V10.10 - Add class variables and methods to be used to correctly calculate the probabilities of single pixel clusters
78 // V10.11 - Allow subdetector ID=5 for FPix R2P2 [allows better internal labeling of templates]
79 // V10.12 - Enforce minimum signal size in pixel charge uncertainty calculation
80 // V10.13 - Update the variable size [SI_PIXEL_TEMPLATE_USE_BOOST] option so that it works with VI's enhancements
81 // V10.20 - Add directory path selection to the ascii pushfile method
82 // V10.21 - Address runtime issues in pushfile() for gcc 7.X due to using tempfile as char string + misc. cleanup [Petar]
83 // V10.22 - Move templateStore to the heap, fix variable name in pushfile() [Petar]
84 // V10.24 - Add sideload() + associated gymnastics [Petar and Oz]
85 
86 // Created by Morris Swartz on 10/27/06.
87 //
88 //
89 
90 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
91 #include <cmath>
92 #else
93 #include <math.h>
94 #endif
95 #include <algorithm>
96 #include <vector>
97 #include "boost/multi_array.hpp"
98 #include <iostream>
99 #include <iomanip>
100 #include <sstream>
101 #include <fstream>
102 #include <list>
103 
104 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
105 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplate.h"
106 #include "CondFormats/SiPixelTransient/interface/SimplePixel.h"
107 #include "FWCore/Utilities/interface/FileInPath.h"
108 #include "FWCore/MessageLogger/interface/MessageLogger.h"
109 #include "FWCore/Utilities/interface/isFinite.h"
110 #define LOGERROR(x) LogError(x)
111 #define LOGINFO(x) LogInfo(x)
112 #define LOGWARNING(x) LogWarning(x)
113 #define ENDL " "
114 #include "FWCore/Utilities/interface/Exception.h"
115 using namespace edm;
116 #else
117 #include "SiPixelTemplate.h"
118 #include "SimplePixel.h"
119 #define LOGERROR(x) std::cout << x << ": "
120 #define LOGINFO(x) std::cout << x << ": "
121 #define ENDL std::endl
122 #endif
123 
124 //****************************************************************
129 //****************************************************************
130 bool SiPixelTemplate::pushfile(int filenum, std::vector<SiPixelTemplateStore>& pixelTemp, std::string dir) {
131  // Add template stored in external file numbered filenum to theTemplateStore
132 
133  // Local variables
134  int i, j, k, l;
135  float qavg_avg;
136  char c;
137  const int code_version = {17};
138 
139  // Create a filename for this run
140  std::string tempfile = std::to_string(filenum);
141 
142 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
143  // If integer filenum has less than 4 digits, prepend 0's until we have four numerical characters, e.g. "0292"
144  int nzeros = 4 - tempfile.length();
145  if (nzeros > 0)
146  tempfile = std::string(nzeros, '0') + tempfile;
148 
149  tempfile = dir + "template_summary_zp" + tempfile + ".out";
150  edm::FileInPath file(tempfile); // Find the file in CMSSW
151  tempfile = file.fullPath(); // Put it back with the whole path.
152 
153 #else
154  // This is the same as above, but more elegant. (Elegance not allowed in CMSSW...)
155  std::ostringstream tout;
156  tout << "template_summary_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
157  tempfile = tout.str();
158 
159 #endif
160 
161  // Open the template file
162  //
163  std::ifstream in_file(tempfile);
164  if (in_file.is_open() && in_file.good()) {
165  // Create a local template storage entry
166 
167  SiPixelTemplateStore theCurrentTemp;
168 
169  // Read-in a header string first and print it
170 
171  for (i = 0; (c = in_file.get()) != '\n'; ++i) {
172  if (i < 79) {
173  theCurrentTemp.head.title[i] = c;
174  }
175  }
176  if (i > 78) {
177  i = 78;
178  }
179  theCurrentTemp.head.title[i + 1] = '\0';
180  LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
181 
182  // next, the header information
183 
184  in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
185  theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
186  theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
187  theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
188  theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
189  theCurrentTemp.head.zsize;
190 
191  if (in_file.fail()) {
192  LOGERROR("SiPixelTemplate") << "Error reading file 0A, no template load" << ENDL;
193  return false;
194  }
195 
196  if (theCurrentTemp.head.templ_version > 17) {
197  in_file >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
198  theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
199 
200  if (in_file.fail()) {
201  LOGERROR("SiPixelTemplate") << "Error reading file 0B, no template load" << ENDL;
202  return false;
203  }
204  } else {
205  theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
206  theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth / 2.f;
207  theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth / 2.f;
208  theCurrentTemp.head.fbin[0] = 1.5f;
209  theCurrentTemp.head.fbin[1] = 1.00f;
210  theCurrentTemp.head.fbin[2] = 0.85f;
211  }
212 
213  LOGINFO("SiPixelTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
214  << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
215  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
216  << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
217  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
218  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
219  << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
220  << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold " << theCurrentTemp.head.ss50
221  << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", y Lorentz Bias "
222  << theCurrentTemp.head.lorybias << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
223  << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
224  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
225  << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
226  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
227  << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
228 
229  if (theCurrentTemp.head.templ_version < code_version) {
230  LOGERROR("SiPixelTemplate") << "code expects version " << code_version << " finds "
231  << theCurrentTemp.head.templ_version << ", no template load" << ENDL;
232  return false;
233  }
234 
235 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
236 
237  // next, layout the 1-d/2-d structures needed to store template
238 
239  theCurrentTemp.cotbetaY = std::vector<float>(theCurrentTemp.head.NTy);
240  theCurrentTemp.cotbetaX = std::vector<float>(theCurrentTemp.head.NTyx);
241  theCurrentTemp.cotalphaX = std::vector<float>(theCurrentTemp.head.NTxx);
242 
243  theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
244 
245  theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
246 
247 #endif
248 
249  // next, loop over all y-angle entries
250 
251  for (i = 0; i < theCurrentTemp.head.NTy; ++i) {
252  in_file >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0] >>
253  theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2];
254 
255  if (in_file.fail()) {
256  LOGERROR("SiPixelTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum
257  << ENDL;
258  return false;
259  }
260 
261  // Calculate the alpha, beta, and cot(beta) for this entry
262 
263  theCurrentTemp.enty[i].alpha =
264  static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
265 
266  theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0] / theCurrentTemp.enty[i].costrk[2];
267 
268  theCurrentTemp.enty[i].beta =
269  static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
270 
271  theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1] / theCurrentTemp.enty[i].costrk[2];
272 
273  in_file >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].symax >>
274  theCurrentTemp.enty[i].dyone >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].sxmax >>
275  theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
276 
277  if (in_file.fail()) {
278  LOGERROR("SiPixelTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum
279  << ENDL;
280  return false;
281  }
282 
283  in_file >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo >>
284  theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clsleny >>
285  theCurrentTemp.enty[i].clslenx;
286 
287  if (in_file.fail()) {
288  LOGERROR("SiPixelTemplate") << "Error reading file 3, no template load, run # " << theCurrentTemp.enty[i].runnum
289  << ENDL;
290  return false;
291  }
292 
293  if (theCurrentTemp.enty[i].qmin <= 0.) {
294  LOGERROR("SiPixelTemplate") << "Error in template ID " << theCurrentTemp.head.ID
295  << " qmin = " << theCurrentTemp.enty[i].qmin << ", run # "
296  << theCurrentTemp.enty[i].runnum << ENDL;
297  return false;
298  }
299 
300  for (j = 0; j < 2; ++j) {
301  in_file >> theCurrentTemp.enty[i].ypar[j][0] >> theCurrentTemp.enty[i].ypar[j][1] >>
302  theCurrentTemp.enty[i].ypar[j][2] >> theCurrentTemp.enty[i].ypar[j][3] >> theCurrentTemp.enty[i].ypar[j][4];
303 
304  if (in_file.fail()) {
305  LOGERROR("SiPixelTemplate") << "Error reading file 4, no template load, run # "
306  << theCurrentTemp.enty[i].runnum << ENDL;
307  return false;
308  }
309  }
310 
311  for (j = 0; j < 9; ++j) {
312  for (k = 0; k < TYSIZE; ++k) {
313  in_file >> theCurrentTemp.enty[i].ytemp[j][k];
314  }
315 
316  if (in_file.fail()) {
317  LOGERROR("SiPixelTemplate") << "Error reading file 5, no template load, run # "
318  << theCurrentTemp.enty[i].runnum << ENDL;
319  return false;
320  }
321  }
322 
323  for (j = 0; j < 2; ++j) {
324  in_file >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1] >>
325  theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
326 
327  if (in_file.fail()) {
328  LOGERROR("SiPixelTemplate") << "Error reading file 6, no template load, run # "
329  << theCurrentTemp.enty[i].runnum << ENDL;
330  return false;
331  }
332  }
333 
334  qavg_avg = 0.f;
335  for (j = 0; j < 9; ++j) {
336  for (k = 0; k < TXSIZE; ++k) {
337  in_file >> theCurrentTemp.enty[i].xtemp[j][k];
338  qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];
339  }
340 
341  if (in_file.fail()) {
342  LOGERROR("SiPixelTemplate") << "Error reading file 7, no template load, run # "
343  << theCurrentTemp.enty[i].runnum << ENDL;
344  return false;
345  }
346  }
347  theCurrentTemp.enty[i].qavg_avg = qavg_avg / 9.;
348 
349  for (j = 0; j < 4; ++j) {
350  in_file >> theCurrentTemp.enty[i].yavg[j] >> theCurrentTemp.enty[i].yrms[j] >> theCurrentTemp.enty[i].ygx0[j] >>
351  theCurrentTemp.enty[i].ygsig[j];
352 
353  if (in_file.fail()) {
354  LOGERROR("SiPixelTemplate") << "Error reading file 8, no template load, run # "
355  << theCurrentTemp.enty[i].runnum << ENDL;
356  return false;
357  }
358  }
359 
360  for (j = 0; j < 4; ++j) {
361  in_file >> theCurrentTemp.enty[i].yflpar[j][0] >> theCurrentTemp.enty[i].yflpar[j][1] >>
362  theCurrentTemp.enty[i].yflpar[j][2] >> theCurrentTemp.enty[i].yflpar[j][3] >>
363  theCurrentTemp.enty[i].yflpar[j][4] >> theCurrentTemp.enty[i].yflpar[j][5];
364 
365  if (in_file.fail()) {
366  LOGERROR("SiPixelTemplate") << "Error reading file 9, no template load, run # "
367  << theCurrentTemp.enty[i].runnum << ENDL;
368  return false;
369  }
370  }
371 
372  for (j = 0; j < 4; ++j) {
373  in_file >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >>
374  theCurrentTemp.enty[i].xgsig[j];
375 
376  if (in_file.fail()) {
377  LOGERROR("SiPixelTemplate") << "Error reading file 10, no template load, run # "
378  << theCurrentTemp.enty[i].runnum << ENDL;
379  return false;
380  }
381  }
382 
383  for (j = 0; j < 4; ++j) {
384  in_file >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >>
385  theCurrentTemp.enty[i].xflpar[j][2] >> theCurrentTemp.enty[i].xflpar[j][3] >>
386  theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
387 
388  if (in_file.fail()) {
389  LOGERROR("SiPixelTemplate") << "Error reading file 11, no template load, run # "
390  << theCurrentTemp.enty[i].runnum << ENDL;
391  return false;
392  }
393  }
394 
395  for (j = 0; j < 4; ++j) {
396  in_file >> theCurrentTemp.enty[i].chi2yavg[j] >> theCurrentTemp.enty[i].chi2ymin[j] >>
397  theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j];
398 
399  if (in_file.fail()) {
400  LOGERROR("SiPixelTemplate") << "Error reading file 12, no template load, run # "
401  << theCurrentTemp.enty[i].runnum << ENDL;
402  return false;
403  }
404  }
405 
406  for (j = 0; j < 4; ++j) {
407  in_file >> theCurrentTemp.enty[i].yavgc2m[j] >> theCurrentTemp.enty[i].yrmsc2m[j] >>
408  theCurrentTemp.enty[i].chi2yavgc2m[j] >> theCurrentTemp.enty[i].chi2yminc2m[j];
409 
410  if (in_file.fail()) {
411  LOGERROR("SiPixelTemplate") << "Error reading file 13, no template load, run # "
412  << theCurrentTemp.enty[i].runnum << ENDL;
413  return false;
414  }
415  }
416 
417  for (j = 0; j < 4; ++j) {
418  in_file >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >>
419  theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
420 
421  if (in_file.fail()) {
422  LOGERROR("SiPixelTemplate") << "Error reading file 14, no template load, run # "
423  << theCurrentTemp.enty[i].runnum << ENDL;
424  return false;
425  }
426  }
427 
428  for (j = 0; j < 4; ++j) {
429  in_file >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >>
430  theCurrentTemp.enty[i].ygx0gen[j] >> theCurrentTemp.enty[i].ygsiggen[j];
431 
432  if (in_file.fail()) {
433  LOGERROR("SiPixelTemplate") << "Error reading file 14a, no template load, run # "
434  << theCurrentTemp.enty[i].runnum << ENDL;
435  return false;
436  }
437  }
438 
439  for (j = 0; j < 4; ++j) {
440  in_file >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >>
441  theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
442 
443  if (in_file.fail()) {
444  LOGERROR("SiPixelTemplate") << "Error reading file 14b, no template load, run # "
445  << theCurrentTemp.enty[i].runnum << ENDL;
446  return false;
447  }
448  }
449 
450  in_file >> theCurrentTemp.enty[i].chi2yavgone >> theCurrentTemp.enty[i].chi2yminone >>
451  theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2 >>
452  theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav >>
453  theCurrentTemp.enty[i].r_qMeas_qTrue >> theCurrentTemp.enty[i].spare[0];
454 
455  if (in_file.fail()) {
456  LOGERROR("SiPixelTemplate") << "Error reading file 15, no template load, run # "
457  << theCurrentTemp.enty[i].runnum << ENDL;
458  return false;
459  }
460 
461  in_file >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >>
462  theCurrentTemp.enty[i].kappavav2 >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1] >>
463  theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracyone >> theCurrentTemp.enty[i].fracxone >>
464  theCurrentTemp.enty[i].fracytwo >> theCurrentTemp.enty[i].fracxtwo;
465  // theCurrentTemp.enty[i].qbfrac[3] = 1. - theCurrentTemp.enty[i].qbfrac[0] - theCurrentTemp.enty[i].qbfrac[1] - theCurrentTemp.enty[i].qbfrac[2];
466 
467  if (in_file.fail()) {
468  LOGERROR("SiPixelTemplate") << "Error reading file 16, no template load, run # "
469  << theCurrentTemp.enty[i].runnum << ENDL;
470  return false;
471  }
472  }
473 
474  // next, loop over all barrel x-angle entries
475 
476  for (k = 0; k < theCurrentTemp.head.NTyx; ++k) {
477  for (i = 0; i < theCurrentTemp.head.NTxx; ++i) {
478  in_file >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0] >>
479  theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2];
480 
481  if (in_file.fail()) {
482  LOGERROR("SiPixelTemplate") << "Error reading file 17, no template load, run # "
483  << theCurrentTemp.entx[k][i].runnum << ENDL;
484  return false;
485  }
486 
487  // Calculate the alpha, beta, and cot(beta) for this entry
488 
489  theCurrentTemp.entx[k][i].alpha = static_cast<float>(
490  atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
491 
492  theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0] / theCurrentTemp.entx[k][i].costrk[2];
493 
494  theCurrentTemp.entx[k][i].beta = static_cast<float>(
495  atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
496 
497  theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1] / theCurrentTemp.entx[k][i].costrk[2];
498 
499  in_file >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >>
500  theCurrentTemp.entx[k][i].symax >> theCurrentTemp.entx[k][i].dyone >> theCurrentTemp.entx[k][i].syone >>
501  theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
502 
503  if (in_file.fail()) {
504  LOGERROR("SiPixelTemplate") << "Error reading file 18, no template load, run # "
505  << theCurrentTemp.entx[k][i].runnum << ENDL;
506  return false;
507  }
508 
509  in_file >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >>
510  theCurrentTemp.entx[k][i].dxtwo >> theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >>
511  theCurrentTemp.entx[k][i].clsleny >> theCurrentTemp.entx[k][i].clslenx;
512  // >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
513 
514  if (in_file.fail()) {
515  LOGERROR("SiPixelTemplate") << "Error reading file 19, no template load, run # "
516  << theCurrentTemp.entx[k][i].runnum << ENDL;
517  return false;
518  }
519 
520  for (j = 0; j < 2; ++j) {
521  in_file >> theCurrentTemp.entx[k][i].ypar[j][0] >> theCurrentTemp.entx[k][i].ypar[j][1] >>
522  theCurrentTemp.entx[k][i].ypar[j][2] >> theCurrentTemp.entx[k][i].ypar[j][3] >>
523  theCurrentTemp.entx[k][i].ypar[j][4];
524 
525  if (in_file.fail()) {
526  LOGERROR("SiPixelTemplate") << "Error reading file 20, no template load, run # "
527  << theCurrentTemp.entx[k][i].runnum << ENDL;
528  return false;
529  }
530  }
531 
532  for (j = 0; j < 9; ++j) {
533  for (l = 0; l < TYSIZE; ++l) {
534  in_file >> theCurrentTemp.entx[k][i].ytemp[j][l];
535  }
536 
537  if (in_file.fail()) {
538  LOGERROR("SiPixelTemplate") << "Error reading file 21, no template load, run # "
539  << theCurrentTemp.entx[k][i].runnum << ENDL;
540  return false;
541  }
542  }
543 
544  for (j = 0; j < 2; ++j) {
545  in_file >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1] >>
546  theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >>
547  theCurrentTemp.entx[k][i].xpar[j][4];
548 
549  if (in_file.fail()) {
550  LOGERROR("SiPixelTemplate") << "Error reading file 22, no template load, run # "
551  << theCurrentTemp.entx[k][i].runnum << ENDL;
552  return false;
553  }
554  }
555 
556  qavg_avg = 0.f;
557  for (j = 0; j < 9; ++j) {
558  for (l = 0; l < TXSIZE; ++l) {
559  in_file >> theCurrentTemp.entx[k][i].xtemp[j][l];
560  qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];
561  }
562 
563  if (in_file.fail()) {
564  LOGERROR("SiPixelTemplate") << "Error reading file 23, no template load, run # "
565  << theCurrentTemp.entx[k][i].runnum << ENDL;
566  return false;
567  }
568  }
569  theCurrentTemp.entx[k][i].qavg_avg = qavg_avg / 9.;
570 
571  for (j = 0; j < 4; ++j) {
572  in_file >> theCurrentTemp.entx[k][i].yavg[j] >> theCurrentTemp.entx[k][i].yrms[j] >>
573  theCurrentTemp.entx[k][i].ygx0[j] >> theCurrentTemp.entx[k][i].ygsig[j];
574 
575  if (in_file.fail()) {
576  LOGERROR("SiPixelTemplate") << "Error reading file 24, no template load, run # "
577  << theCurrentTemp.entx[k][i].runnum << ENDL;
578  return false;
579  }
580  }
581 
582  for (j = 0; j < 4; ++j) {
583  in_file >> theCurrentTemp.entx[k][i].yflpar[j][0] >> theCurrentTemp.entx[k][i].yflpar[j][1] >>
584  theCurrentTemp.entx[k][i].yflpar[j][2] >> theCurrentTemp.entx[k][i].yflpar[j][3] >>
585  theCurrentTemp.entx[k][i].yflpar[j][4] >> theCurrentTemp.entx[k][i].yflpar[j][5];
586 
587  if (in_file.fail()) {
588  LOGERROR("SiPixelTemplate") << "Error reading file 25, no template load, run # "
589  << theCurrentTemp.entx[k][i].runnum << ENDL;
590  return false;
591  }
592  }
593 
594  for (j = 0; j < 4; ++j) {
595  in_file >> theCurrentTemp.entx[k][i].xavg[j] >> theCurrentTemp.entx[k][i].xrms[j] >>
596  theCurrentTemp.entx[k][i].xgx0[j] >> theCurrentTemp.entx[k][i].xgsig[j];
597 
598  if (in_file.fail()) {
599  LOGERROR("SiPixelTemplate") << "Error reading file 26, no template load, run # "
600  << theCurrentTemp.entx[k][i].runnum << ENDL;
601  return false;
602  }
603  }
604 
605  for (j = 0; j < 4; ++j) {
606  in_file >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >>
607  theCurrentTemp.entx[k][i].xflpar[j][2] >> theCurrentTemp.entx[k][i].xflpar[j][3] >>
608  theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
609 
610  if (in_file.fail()) {
611  LOGERROR("SiPixelTemplate") << "Error reading file 27, no template load, run # "
612  << theCurrentTemp.entx[k][i].runnum << ENDL;
613  return false;
614  }
615  }
616 
617  for (j = 0; j < 4; ++j) {
618  in_file >> theCurrentTemp.entx[k][i].chi2yavg[j] >> theCurrentTemp.entx[k][i].chi2ymin[j] >>
619  theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j];
620 
621  if (in_file.fail()) {
622  LOGERROR("SiPixelTemplate") << "Error reading file 28, no template load, run # "
623  << theCurrentTemp.entx[k][i].runnum << ENDL;
624  return false;
625  }
626  }
627 
628  for (j = 0; j < 4; ++j) {
629  in_file >> theCurrentTemp.entx[k][i].yavgc2m[j] >> theCurrentTemp.entx[k][i].yrmsc2m[j] >>
630  theCurrentTemp.entx[k][i].chi2yavgc2m[j] >> theCurrentTemp.entx[k][i].chi2yminc2m[j];
631 
632  if (in_file.fail()) {
633  LOGERROR("SiPixelTemplate") << "Error reading file 29, no template load, run # "
634  << theCurrentTemp.entx[k][i].runnum << ENDL;
635  return false;
636  }
637  }
638 
639  for (j = 0; j < 4; ++j) {
640  in_file >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >>
641  theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
642 
643  if (in_file.fail()) {
644  LOGERROR("SiPixelTemplate") << "Error reading file 30, no template load, run # "
645  << theCurrentTemp.entx[k][i].runnum << ENDL;
646  return false;
647  }
648  }
649 
650  for (j = 0; j < 4; ++j) {
651  in_file >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >>
652  theCurrentTemp.entx[k][i].ygx0gen[j] >> theCurrentTemp.entx[k][i].ygsiggen[j];
653 
654  if (in_file.fail()) {
655  LOGERROR("SiPixelTemplate") << "Error reading file 30a, no template load, run # "
656  << theCurrentTemp.entx[k][i].runnum << ENDL;
657  return false;
658  }
659  }
660 
661  for (j = 0; j < 4; ++j) {
662  in_file >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j] >>
663  theCurrentTemp.entx[k][i].xgx0gen[j] >> theCurrentTemp.entx[k][i].xgsiggen[j];
664 
665  if (in_file.fail()) {
666  LOGERROR("SiPixelTemplate") << "Error reading file 30b, no template load, run # "
667  << theCurrentTemp.entx[k][i].runnum << ENDL;
668  return false;
669  }
670  }
671 
672  in_file >> theCurrentTemp.entx[k][i].chi2yavgone >> theCurrentTemp.entx[k][i].chi2yminone >>
673  theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >>
674  theCurrentTemp.entx[k][i].qmin2 >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >>
675  theCurrentTemp.entx[k][i].kappavav >> theCurrentTemp.entx[k][i].r_qMeas_qTrue >>
676  theCurrentTemp.entx[k][i].spare[0];
677 
678  if (in_file.fail()) {
679  LOGERROR("SiPixelTemplate") << "Error reading file 31, no template load, run # "
680  << theCurrentTemp.entx[k][i].runnum << ENDL;
681  return false;
682  }
683 
684  in_file >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >>
685  theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].qbfrac[0] >>
686  theCurrentTemp.entx[k][i].qbfrac[1] >> theCurrentTemp.entx[k][i].qbfrac[2] >>
687  theCurrentTemp.entx[k][i].fracyone >> theCurrentTemp.entx[k][i].fracxone >>
688  theCurrentTemp.entx[k][i].fracytwo >> theCurrentTemp.entx[k][i].fracxtwo;
689  // theCurrentTemp.entx[k][i].qbfrac[3] = 1. - theCurrentTemp.entx[k][i].qbfrac[0] - theCurrentTemp.entx[k][i].qbfrac[1] - theCurrentTemp.entx[k][i].qbfrac[2];
690 
691  if (in_file.fail()) {
692  LOGERROR("SiPixelTemplate") << "Error reading file 32, no template load, run # "
693  << theCurrentTemp.entx[k][i].runnum << ENDL;
694  return false;
695  }
696  }
697  }
698 
699  in_file.close();
700 
701  // Add this template to the store
702 
703  pixelTemp.push_back(theCurrentTemp);
704 
705  postInit(pixelTemp);
706  return true;
707 
708  } else {
709  // If file didn't open, report this
710 
711  LOGERROR("SiPixelTemplate") << "Error opening File" << tempfile << ENDL;
712  return false;
713  }
714 
715 } // TempInit
716 
717 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
718 
719 //****************************************************************
723 //****************************************************************
724 bool SiPixelTemplate::pushfile(const SiPixelTemplateDBObject& dbobject, std::vector<SiPixelTemplateStore>& pixelTemp) {
725  // Add template stored in external dbobject to theTemplateStore
726 
727  // Local variables
728  int i, j, k, l;
729  float qavg_avg;
730  const int code_version = {17};
731 
732  // We must create a new object because dbobject must be a const and our stream must not be
733  auto db(dbobject.reader());
734 
735  // Create a local template storage entry
737  auto tmpPtr = std::make_unique<SiPixelTemplateStore>(); // must be allocated on the heap instead
738  auto& theCurrentTemp = *tmpPtr;
739 
740  // Fill the template storage for each template calibration stored in the db
741  for (int m = 0; m < db.numOfTempl(); ++m) {
742  // Read-in a header string first and print it
743 
744  SiPixelTemplateDBObject::char2float temp;
745  for (i = 0; i < 20; ++i) {
746  temp.f = db.sVector()[db.index()];
747  theCurrentTemp.head.title[4 * i] = temp.c[0];
748  theCurrentTemp.head.title[4 * i + 1] = temp.c[1];
749  theCurrentTemp.head.title[4 * i + 2] = temp.c[2];
750  theCurrentTemp.head.title[4 * i + 3] = temp.c[3];
751  db.incrementIndex(1);
752  }
753  theCurrentTemp.head.title[79] = '\0';
754  LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
755 
756  // next, the header information
757 
758  db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
759  theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
760  theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
761  theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
762  theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
763  theCurrentTemp.head.zsize;
764 
765  if (db.fail()) {
766  LOGERROR("SiPixelTemplate") << "Error reading file 0A, no template load" << ENDL;
767  return false;
768  }
769 
770  LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title
771  << " code version = " << code_version << " object version "
772  << theCurrentTemp.head.templ_version << ENDL;
773 
774  if (theCurrentTemp.head.templ_version > 17) {
775  db >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
776  theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
777 
778  if (db.fail()) {
779  LOGERROR("SiPixelTemplate") << "Error reading file 0B, no template load" << ENDL;
780  return false;
781  }
782  } else {
783  theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
784  theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth / 2.f;
785  theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth / 2.f;
786  theCurrentTemp.head.fbin[0] = 1.50f;
787  theCurrentTemp.head.fbin[1] = 1.00f;
788  theCurrentTemp.head.fbin[2] = 0.85f;
789  //std::cout<<" set fbin "<< theCurrentTemp.head.fbin[0]<<" "<<theCurrentTemp.head.fbin[1]<<" "
790  // <<theCurrentTemp.head.fbin[2]<<std::endl;
791  }
792 
793  LOGINFO("SiPixelTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
794  << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
795  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
796  << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
797  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
798  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
799  << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
800  << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold " << theCurrentTemp.head.ss50
801  << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", y Lorentz Bias "
802  << theCurrentTemp.head.lorybias << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
803  << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
804  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
805  << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
806  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
807  << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
808 
809  if (theCurrentTemp.head.templ_version < code_version) {
810  LOGINFO("SiPixelTemplate") << "code expects version " << code_version << " finds "
811  << theCurrentTemp.head.templ_version << ", load anyway " << ENDL;
812  //return false; // dk
813  }
814 
815 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
816 
817  // next, layout the 1-d/2-d structures needed to store template
818  theCurrentTemp.cotbetaY = std::vector<float>(theCurrentTemp.head.NTy);
819  theCurrentTemp.cotbetaX = std::vector<float>(theCurrentTemp.head.NTyx);
820  theCurrentTemp.cotalphaX = std::vector<float>(theCurrentTemp.head.NTxx);
821  theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
822  theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
823 
824 #endif
825 
826  // next, loop over all barrel y-angle entries
827 
828  for (i = 0; i < theCurrentTemp.head.NTy; ++i) {
829  db >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0] >> theCurrentTemp.enty[i].costrk[1] >>
830  theCurrentTemp.enty[i].costrk[2];
831 
832  if (db.fail()) {
833  LOGERROR("SiPixelTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum
834  << ENDL;
835  return false;
836  }
837 
838  // Calculate the alpha, beta, and cot(beta) for this entry
839 
840  theCurrentTemp.enty[i].alpha =
841  static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
842 
843  theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0] / theCurrentTemp.enty[i].costrk[2];
844 
845  theCurrentTemp.enty[i].beta =
846  static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
847 
848  theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1] / theCurrentTemp.enty[i].costrk[2];
849 
850  db >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].symax >>
851  theCurrentTemp.enty[i].dyone >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].sxmax >>
852  theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
853 
854  if (db.fail()) {
855  LOGERROR("SiPixelTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum
856  << ENDL;
857  return false;
858  }
859 
860  db >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo >>
861  theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clsleny >>
862  theCurrentTemp.enty[i].clslenx;
863  // >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav;
864 
865  if (db.fail()) {
866  LOGERROR("SiPixelTemplate") << "Error reading file 3, no template load, run # " << theCurrentTemp.enty[i].runnum
867  << ENDL;
868  return false;
869  }
870 
871  if (theCurrentTemp.enty[i].qmin <= 0.) {
872  LOGERROR("SiPixelTemplate") << "Error in template ID " << theCurrentTemp.head.ID
873  << " qmin = " << theCurrentTemp.enty[i].qmin << ", run # "
874  << theCurrentTemp.enty[i].runnum << ENDL;
875  return false;
876  }
877 
878  for (j = 0; j < 2; ++j) {
879  db >> theCurrentTemp.enty[i].ypar[j][0] >> theCurrentTemp.enty[i].ypar[j][1] >>
880  theCurrentTemp.enty[i].ypar[j][2] >> theCurrentTemp.enty[i].ypar[j][3] >> theCurrentTemp.enty[i].ypar[j][4];
881 
882  if (db.fail()) {
883  LOGERROR("SiPixelTemplate") << "Error reading file 4, no template load, run # "
884  << theCurrentTemp.enty[i].runnum << ENDL;
885  return false;
886  }
887  }
888 
889  for (j = 0; j < 9; ++j) {
890  for (k = 0; k < TYSIZE; ++k) {
891  db >> theCurrentTemp.enty[i].ytemp[j][k];
892  }
893 
894  if (db.fail()) {
895  LOGERROR("SiPixelTemplate") << "Error reading file 5, no template load, run # "
896  << theCurrentTemp.enty[i].runnum << ENDL;
897  return false;
898  }
899  }
900 
901  for (j = 0; j < 2; ++j) {
902  db >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1] >>
903  theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
904 
905  if (db.fail()) {
906  LOGERROR("SiPixelTemplate") << "Error reading file 6, no template load, run # "
907  << theCurrentTemp.enty[i].runnum << ENDL;
908  return false;
909  }
910  }
911 
912  qavg_avg = 0.f;
913  for (j = 0; j < 9; ++j) {
914  for (k = 0; k < TXSIZE; ++k) {
915  db >> theCurrentTemp.enty[i].xtemp[j][k];
916  qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];
917  }
918 
919  if (db.fail()) {
920  LOGERROR("SiPixelTemplate") << "Error reading file 7, no template load, run # "
921  << theCurrentTemp.enty[i].runnum << ENDL;
922  return false;
923  }
924  }
925  theCurrentTemp.enty[i].qavg_avg = qavg_avg / 9.;
926 
927  for (j = 0; j < 4; ++j) {
928  db >> theCurrentTemp.enty[i].yavg[j] >> theCurrentTemp.enty[i].yrms[j] >> theCurrentTemp.enty[i].ygx0[j] >>
929  theCurrentTemp.enty[i].ygsig[j];
930 
931  if (db.fail()) {
932  LOGERROR("SiPixelTemplate") << "Error reading file 8, no template load, run # "
933  << theCurrentTemp.enty[i].runnum << ENDL;
934  return false;
935  }
936  }
937 
938  for (j = 0; j < 4; ++j) {
939  db >> theCurrentTemp.enty[i].yflpar[j][0] >> theCurrentTemp.enty[i].yflpar[j][1] >>
940  theCurrentTemp.enty[i].yflpar[j][2] >> theCurrentTemp.enty[i].yflpar[j][3] >>
941  theCurrentTemp.enty[i].yflpar[j][4] >> theCurrentTemp.enty[i].yflpar[j][5];
942 
943  if (db.fail()) {
944  LOGERROR("SiPixelTemplate") << "Error reading file 9, no template load, run # "
945  << theCurrentTemp.enty[i].runnum << ENDL;
946  return false;
947  }
948  }
949 
950  for (j = 0; j < 4; ++j) {
951  db >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >>
952  theCurrentTemp.enty[i].xgsig[j];
953 
954  if (db.fail()) {
955  LOGERROR("SiPixelTemplate") << "Error reading file 10, no template load, run # "
956  << theCurrentTemp.enty[i].runnum << ENDL;
957  return false;
958  }
959  }
960 
961  for (j = 0; j < 4; ++j) {
962  db >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >>
963  theCurrentTemp.enty[i].xflpar[j][2] >> theCurrentTemp.enty[i].xflpar[j][3] >>
964  theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
965 
966  if (db.fail()) {
967  LOGERROR("SiPixelTemplate") << "Error reading file 11, no template load, run # "
968  << theCurrentTemp.enty[i].runnum << ENDL;
969  return false;
970  }
971  }
972 
973  for (j = 0; j < 4; ++j) {
974  db >> theCurrentTemp.enty[i].chi2yavg[j] >> theCurrentTemp.enty[i].chi2ymin[j] >>
975  theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j];
976 
977  if (db.fail()) {
978  LOGERROR("SiPixelTemplate") << "Error reading file 12, no template load, run # "
979  << theCurrentTemp.enty[i].runnum << ENDL;
980  return false;
981  }
982  }
983 
984  for (j = 0; j < 4; ++j) {
985  db >> theCurrentTemp.enty[i].yavgc2m[j] >> theCurrentTemp.enty[i].yrmsc2m[j] >>
986  theCurrentTemp.enty[i].chi2yavgc2m[j] >> theCurrentTemp.enty[i].chi2yminc2m[j];
987 
988  if (db.fail()) {
989  LOGERROR("SiPixelTemplate") << "Error reading file 13, no template load, run # "
990  << theCurrentTemp.enty[i].runnum << ENDL;
991  return false;
992  }
993  }
994 
995  for (j = 0; j < 4; ++j) {
996  db >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >>
997  theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
998 
999  if (db.fail()) {
1000  LOGERROR("SiPixelTemplate") << "Error reading file 14, no template load, run # "
1001  << theCurrentTemp.enty[i].runnum << ENDL;
1002  return false;
1003  }
1004  }
1005 
1006  for (j = 0; j < 4; ++j) {
1007  db >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >>
1008  theCurrentTemp.enty[i].ygx0gen[j] >> theCurrentTemp.enty[i].ygsiggen[j];
1009 
1010  if (db.fail()) {
1011  LOGERROR("SiPixelTemplate") << "Error reading file 14a, no template load, run # "
1012  << theCurrentTemp.enty[i].runnum << ENDL;
1013  return false;
1014  }
1015  }
1016 
1017  for (j = 0; j < 4; ++j) {
1018  db >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >>
1019  theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
1020 
1021  if (db.fail()) {
1022  LOGERROR("SiPixelTemplate") << "Error reading file 14b, no template load, run # "
1023  << theCurrentTemp.enty[i].runnum << ENDL;
1024  return false;
1025  }
1026  }
1027 
1028  db >> theCurrentTemp.enty[i].chi2yavgone >> theCurrentTemp.enty[i].chi2yminone >>
1029  theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2 >>
1030  theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav >>
1031  theCurrentTemp.enty[i].r_qMeas_qTrue >> theCurrentTemp.enty[i].spare[0];
1032 
1033  if (db.fail()) {
1034  LOGERROR("SiPixelTemplate") << "Error reading file 15, no template load, run # "
1035  << theCurrentTemp.enty[i].runnum << ENDL;
1036  return false;
1037  }
1038 
1039  db >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >> theCurrentTemp.enty[i].kappavav2 >>
1040  theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1] >> theCurrentTemp.enty[i].qbfrac[2] >>
1041  theCurrentTemp.enty[i].fracyone >> theCurrentTemp.enty[i].fracxone >> theCurrentTemp.enty[i].fracytwo >>
1042  theCurrentTemp.enty[i].fracxtwo;
1043  // theCurrentTemp.enty[i].qbfrac[3] = 1. - theCurrentTemp.enty[i].qbfrac[0] - theCurrentTemp.enty[i].qbfrac[1] - theCurrentTemp.enty[i].qbfrac[2];
1044 
1045  if (db.fail()) {
1046  LOGERROR("SiPixelTemplate") << "Error reading file 16, no template load, run # "
1047  << theCurrentTemp.enty[i].runnum << ENDL;
1048  return false;
1049  }
1050  }
1051 
1052  // next, loop over all barrel x-angle entries
1053 
1054  for (k = 0; k < theCurrentTemp.head.NTyx; ++k) {
1055  for (i = 0; i < theCurrentTemp.head.NTxx; ++i) {
1056  db >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0] >>
1057  theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2];
1058 
1059  if (db.fail()) {
1060  LOGERROR("SiPixelTemplate") << "Error reading file 17, no template load, run # "
1061  << theCurrentTemp.entx[k][i].runnum << ENDL;
1062  return false;
1063  }
1064 
1065  // Calculate the alpha, beta, and cot(beta) for this entry
1066 
1067  theCurrentTemp.entx[k][i].alpha = static_cast<float>(
1068  atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
1069 
1070  theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0] / theCurrentTemp.entx[k][i].costrk[2];
1071 
1072  theCurrentTemp.entx[k][i].beta = static_cast<float>(
1073  atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
1074 
1075  theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1] / theCurrentTemp.entx[k][i].costrk[2];
1076 
1077  db >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >> theCurrentTemp.entx[k][i].symax >>
1078  theCurrentTemp.entx[k][i].dyone >> theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].sxmax >>
1079  theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
1080 
1081  if (db.fail()) {
1082  LOGERROR("SiPixelTemplate") << "Error reading file 18, no template load, run # "
1083  << theCurrentTemp.entx[k][i].runnum << ENDL;
1084  return false;
1085  }
1086 
1087  db >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >> theCurrentTemp.entx[k][i].dxtwo >>
1088  theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clsleny >>
1089  theCurrentTemp.entx[k][i].clslenx;
1090  // >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
1091 
1092  if (db.fail()) {
1093  LOGERROR("SiPixelTemplate") << "Error reading file 19, no template load, run # "
1094  << theCurrentTemp.entx[k][i].runnum << ENDL;
1095  return false;
1096  }
1097 
1098  for (j = 0; j < 2; ++j) {
1099  db >> theCurrentTemp.entx[k][i].ypar[j][0] >> theCurrentTemp.entx[k][i].ypar[j][1] >>
1100  theCurrentTemp.entx[k][i].ypar[j][2] >> theCurrentTemp.entx[k][i].ypar[j][3] >>
1101  theCurrentTemp.entx[k][i].ypar[j][4];
1102 
1103  if (db.fail()) {
1104  LOGERROR("SiPixelTemplate") << "Error reading file 20, no template load, run # "
1105  << theCurrentTemp.entx[k][i].runnum << ENDL;
1106  return false;
1107  }
1108  }
1109 
1110  for (j = 0; j < 9; ++j) {
1111  for (l = 0; l < TYSIZE; ++l) {
1112  db >> theCurrentTemp.entx[k][i].ytemp[j][l];
1113  }
1114 
1115  if (db.fail()) {
1116  LOGERROR("SiPixelTemplate") << "Error reading file 21, no template load, run # "
1117  << theCurrentTemp.entx[k][i].runnum << ENDL;
1118  return false;
1119  }
1120  }
1121 
1122  for (j = 0; j < 2; ++j) {
1123  db >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1] >>
1124  theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >>
1125  theCurrentTemp.entx[k][i].xpar[j][4];
1126 
1127  if (db.fail()) {
1128  LOGERROR("SiPixelTemplate") << "Error reading file 22, no template load, run # "
1129  << theCurrentTemp.entx[k][i].runnum << ENDL;
1130  return false;
1131  }
1132  }
1133 
1134  qavg_avg = 0.f;
1135  for (j = 0; j < 9; ++j) {
1136  for (l = 0; l < TXSIZE; ++l) {
1137  db >> theCurrentTemp.entx[k][i].xtemp[j][l];
1138  qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];
1139  }
1140 
1141  if (db.fail()) {
1142  LOGERROR("SiPixelTemplate") << "Error reading file 23, no template load, run # "
1143  << theCurrentTemp.entx[k][i].runnum << ENDL;
1144  return false;
1145  }
1146  }
1147  theCurrentTemp.entx[k][i].qavg_avg = qavg_avg / 9.;
1148 
1149  for (j = 0; j < 4; ++j) {
1150  db >> theCurrentTemp.entx[k][i].yavg[j] >> theCurrentTemp.entx[k][i].yrms[j] >>
1151  theCurrentTemp.entx[k][i].ygx0[j] >> theCurrentTemp.entx[k][i].ygsig[j];
1152 
1153  if (db.fail()) {
1154  LOGERROR("SiPixelTemplate") << "Error reading file 24, no template load, run # "
1155  << theCurrentTemp.entx[k][i].runnum << ENDL;
1156  return false;
1157  }
1158  }
1159 
1160  for (j = 0; j < 4; ++j) {
1161  db >> theCurrentTemp.entx[k][i].yflpar[j][0] >> theCurrentTemp.entx[k][i].yflpar[j][1] >>
1162  theCurrentTemp.entx[k][i].yflpar[j][2] >> theCurrentTemp.entx[k][i].yflpar[j][3] >>
1163  theCurrentTemp.entx[k][i].yflpar[j][4] >> theCurrentTemp.entx[k][i].yflpar[j][5];
1164 
1165  if (db.fail()) {
1166  LOGERROR("SiPixelTemplate") << "Error reading file 25, no template load, run # "
1167  << theCurrentTemp.entx[k][i].runnum << ENDL;
1168  return false;
1169  }
1170  }
1171 
1172  for (j = 0; j < 4; ++j) {
1173  db >> theCurrentTemp.entx[k][i].xavg[j] >> theCurrentTemp.entx[k][i].xrms[j] >>
1174  theCurrentTemp.entx[k][i].xgx0[j] >> theCurrentTemp.entx[k][i].xgsig[j];
1175 
1176  if (db.fail()) {
1177  LOGERROR("SiPixelTemplate") << "Error reading file 26, no template load, run # "
1178  << theCurrentTemp.entx[k][i].runnum << ENDL;
1179  return false;
1180  }
1181  }
1182 
1183  for (j = 0; j < 4; ++j) {
1184  db >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >>
1185  theCurrentTemp.entx[k][i].xflpar[j][2] >> theCurrentTemp.entx[k][i].xflpar[j][3] >>
1186  theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
1187 
1188  if (db.fail()) {
1189  LOGERROR("SiPixelTemplate") << "Error reading file 27, no template load, run # "
1190  << theCurrentTemp.entx[k][i].runnum << ENDL;
1191  return false;
1192  }
1193  }
1194 
1195  for (j = 0; j < 4; ++j) {
1196  db >> theCurrentTemp.entx[k][i].chi2yavg[j] >> theCurrentTemp.entx[k][i].chi2ymin[j] >>
1197  theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j];
1198 
1199  if (db.fail()) {
1200  LOGERROR("SiPixelTemplate") << "Error reading file 28, no template load, run # "
1201  << theCurrentTemp.entx[k][i].runnum << ENDL;
1202  return false;
1203  }
1204  }
1205 
1206  for (j = 0; j < 4; ++j) {
1207  db >> theCurrentTemp.entx[k][i].yavgc2m[j] >> theCurrentTemp.entx[k][i].yrmsc2m[j] >>
1208  theCurrentTemp.entx[k][i].chi2yavgc2m[j] >> theCurrentTemp.entx[k][i].chi2yminc2m[j];
1209 
1210  if (db.fail()) {
1211  LOGERROR("SiPixelTemplate") << "Error reading file 29, no template load, run # "
1212  << theCurrentTemp.entx[k][i].runnum << ENDL;
1213  return false;
1214  }
1215  }
1216 
1217  for (j = 0; j < 4; ++j) {
1218  db >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >>
1219  theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
1220 
1221  if (db.fail()) {
1222  LOGERROR("SiPixelTemplate") << "Error reading file 30, no template load, run # "
1223  << theCurrentTemp.entx[k][i].runnum << ENDL;
1224  return false;
1225  }
1226  }
1227 
1228  for (j = 0; j < 4; ++j) {
1229  db >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >>
1230  theCurrentTemp.entx[k][i].ygx0gen[j] >> theCurrentTemp.entx[k][i].ygsiggen[j];
1231 
1232  if (db.fail()) {
1233  LOGERROR("SiPixelTemplate") << "Error reading file 30a, no template load, run # "
1234  << theCurrentTemp.entx[k][i].runnum << ENDL;
1235  return false;
1236  }
1237  }
1238 
1239  for (j = 0; j < 4; ++j) {
1240  db >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j] >>
1241  theCurrentTemp.entx[k][i].xgx0gen[j] >> theCurrentTemp.entx[k][i].xgsiggen[j];
1242 
1243  if (db.fail()) {
1244  LOGERROR("SiPixelTemplate") << "Error reading file 30b, no template load, run # "
1245  << theCurrentTemp.entx[k][i].runnum << ENDL;
1246  return false;
1247  }
1248  }
1249 
1250  db >> theCurrentTemp.entx[k][i].chi2yavgone >> theCurrentTemp.entx[k][i].chi2yminone >>
1251  theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >>
1252  theCurrentTemp.entx[k][i].qmin2 >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >>
1253  theCurrentTemp.entx[k][i].kappavav >> theCurrentTemp.entx[k][i].r_qMeas_qTrue >>
1254  theCurrentTemp.entx[k][i].spare[0];
1255 
1256  if (db.fail()) {
1257  LOGERROR("SiPixelTemplate") << "Error reading file 31, no template load, run # "
1258  << theCurrentTemp.entx[k][i].runnum << ENDL;
1259  return false;
1260  }
1261 
1262  db >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >>
1263  theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].qbfrac[0] >>
1264  theCurrentTemp.entx[k][i].qbfrac[1] >> theCurrentTemp.entx[k][i].qbfrac[2] >>
1265  theCurrentTemp.entx[k][i].fracyone >> theCurrentTemp.entx[k][i].fracxone >>
1266  theCurrentTemp.entx[k][i].fracytwo >> theCurrentTemp.entx[k][i].fracxtwo;
1267  // theCurrentTemp.entx[k][i].qbfrac[3] = 1. - theCurrentTemp.entx[k][i].qbfrac[0] - theCurrentTemp.entx[k][i].qbfrac[1] - theCurrentTemp.entx[k][i].qbfrac[2];
1268 
1269  if (db.fail()) {
1270  LOGERROR("SiPixelTemplate") << "Error reading file 32, no template load, run # "
1271  << theCurrentTemp.entx[k][i].runnum << ENDL;
1272  return false;
1273  }
1274  }
1275  }
1276 
1277  // Add this template to the store
1278 
1279  pixelTemp.push_back(theCurrentTemp);
1280  }
1281  postInit(pixelTemp);
1282  return true;
1283 
1284 } // TempInit
1285 
1286 #endif
1287 
1288 void SiPixelTemplate::postInit(std::vector<SiPixelTemplateStore>& thePixelTemp_) {
1289  /*
1290  std::cout << "SiPixelTemplate size " << thePixelTemp_.size() << std::endl;
1291  #ifndef SI_PIXEL_TEMPLATE_USE_BOOST
1292  std::cout <<"uses C arrays" << std::endl;
1293  #endif
1294 
1295  int i=0;
1296  for (auto & templ : thePixelTemp_) {
1297  std::cout << i <<':' << templ.head.ID << ' ' << templ.head.NTy <<','<< templ.head.NTyx <<','<< templ.head.NTxx << std::endl;
1298  for ( auto iy=1; iy<templ.head.NTy; ++iy ) { auto & ent = templ.enty[iy]; std::cout << ent.cotbeta <<',' << ent.cotbeta-templ.enty[iy-1].cotbeta << ' '; }
1299  std::cout << std::endl;
1300  for ( auto ix=1; ix<templ.head.NTxx; ++ix ){ auto & ent = templ.entx[0][ix]; std::cout << ent.cotalpha <<','<< ent.cotalpha-templ.entx[0][ix-1].cotalpha << ' ';}
1301  std::cout << std::endl;
1302  ++i;
1303  }
1304  */
1305 
1306  for (auto& templ : thePixelTemp_) {
1307  for (auto iy = 0; iy < templ.head.NTy; ++iy)
1308  templ.cotbetaY[iy] = templ.enty[iy].cotbeta;
1309  for (auto iy = 0; iy < templ.head.NTyx; ++iy)
1310  templ.cotbetaX[iy] = templ.entx[iy][0].cotbeta;
1311  for (auto ix = 0; ix < templ.head.NTxx; ++ix)
1312  templ.cotalphaX[ix] = templ.entx[0][ix].cotalpha;
1313  }
1314 }
1315 
1316 // ************************************************************************************************************
1327 // ************************************************************************************************************
1328 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx) {
1329  // Interpolate for a new set of track angles
1330 
1331  // Local variables
1332  int i, j;
1333  int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx;
1334  float yxratio, xxratio, sxmax, qcorrect, qxtempcor, symax, chi2xavgone, chi2xminone, cota, cotb, cotalpha0, cotbeta0;
1335  float chi2xavg[4], chi2xmin[4], chi2xavgc2m[4], chi2xminc2m[4];
1336 
1337  // If the sideloading is turned on, xtemp_ and ytemp_ are already set to what they need to be.
1338  // So we'll just exit.
1339  if (entry_sideloaded_ != nullptr) {
1340  success_ = true;
1341  return success_;
1342  }
1343 
1344  // Check to see if interpolation is valid
1345  if (id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
1346  cota_current_ = cotalpha;
1347  cotb_current_ = cotbeta;
1348  success_ = true;
1349 
1350  if (id != id_current_) {
1351  // Find the index corresponding to id
1352 
1353  index_id_ = -1;
1354  for (i = 0; i < (int)thePixelTemp_.size(); ++i) {
1355  //std::cout<<i<<" "<<id<<" "<<thePixelTemp_[i].head.ID<<std::endl;
1356 
1357  if (id == thePixelTemp_[i].head.ID) {
1358  index_id_ = i;
1359  id_current_ = id;
1360 
1361  // Copy the detector type to the private variable
1362 
1363  dtype_ = thePixelTemp_[index_id_].head.Dtype;
1364 
1365  // Copy the charge scaling factor to the private variable
1366 
1367  qscale_ = thePixelTemp_[index_id_].head.qscale;
1368 
1369  // Copy the pseudopixel signal size to the private variable
1370 
1371  s50_ = thePixelTemp_[index_id_].head.s50;
1372 
1373  // Copy Qbinning info to private variables
1374 
1375  for (j = 0; j < 3; ++j) {
1376  fbin_[j] = thePixelTemp_[index_id_].head.fbin[j];
1377  }
1378  //std::cout<<" set fbin "<< fbin_[0]<<" "<<fbin_[1]<<" "<<fbin_[2]<<std::endl;
1379 
1380  // Pixel sizes to the private variables
1381 
1382  xsize_ = thePixelTemp_[index_id_].head.xsize;
1383  ysize_ = thePixelTemp_[index_id_].head.ysize;
1384  zsize_ = thePixelTemp_[index_id_].head.zsize;
1385 
1386  break;
1387  }
1388  }
1389  }
1390 
1391 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1392  if (index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
1393  throw cms::Exception("DataCorrupt")
1394  << "SiPixelTemplate::interpolate can't find needed template ID = " << id << std::endl;
1395  }
1396 
1397  //check for nan's
1398  if (!edm::isFinite(cotalpha) || !edm::isFinite(cotbeta)) {
1399  success_ = false;
1400  return success_;
1401  }
1402 #else
1403  assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
1404 #endif
1405 
1406  // qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
1407 
1408  cotalpha0 = thePixelTemp_[index_id_].enty[0].cotalpha;
1409  qcorrect =
1410  std::sqrt((1.f + cotbeta * cotbeta + cotalpha * cotalpha) / (1.f + cotbeta * cotbeta + cotalpha0 * cotalpha0));
1411 
1412  // for some cosmics, the ususal gymnastics are incorrect
1413  cota = cotalpha;
1414  cotb = abs_cotb_ = std::abs(cotbeta);
1415  flip_x_ = false;
1416  flip_y_ = false;
1417  switch (dtype_) {
1418  case 0:
1419  if (cotbeta < 0.f) {
1420  flip_y_ = true;
1421  }
1422  break;
1423  case 1:
1424  if (locBz < 0.f) {
1425  cotb = cotbeta;
1426  } else {
1427  cotb = -cotbeta;
1428  flip_y_ = true;
1429  }
1430  break;
1431  case 2:
1432  case 3:
1433  case 4:
1434  case 5:
1435  if (locBx * locBz < 0.f) {
1436  cota = -cotalpha;
1437  flip_x_ = true;
1438  }
1439  if (locBx > 0.f) {
1440  cotb = cotbeta;
1441  } else {
1442  cotb = -cotbeta;
1443  flip_y_ = true;
1444  }
1445  break;
1446  default:
1447 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1448  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::illegal subdetector ID = " << dtype_ << std::endl;
1449 #else
1450  std::cout << "SiPixelTemplate::illegal subdetector ID = " << dtype_ << std::endl;
1451 #endif
1452  }
1453 
1454  Ny = thePixelTemp_[index_id_].head.NTy;
1455  Nyx = thePixelTemp_[index_id_].head.NTyx;
1456  Nxx = thePixelTemp_[index_id_].head.NTxx;
1457 
1458 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1459  if (Ny < 2 || Nyx < 1 || Nxx < 2) {
1460  throw cms::Exception("DataCorrupt")
1461  << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx
1462  << std::endl;
1463  }
1464 #else
1465  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
1466 #endif
1467  imaxx = Nyx - 1;
1468  imidy = Nxx / 2;
1469 
1470  // next, loop over all y-angle entries
1471 
1472  ilow = 0;
1473  yratio_ = 0.f;
1474 
1475  if (cotb >= thePixelTemp_[index_id_].enty[Ny - 1].cotbeta) {
1476  ilow = Ny - 2;
1477  yratio_ = 1.;
1478  success_ = false;
1479 
1480  } else {
1481  if (cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
1482  for (i = 0; i < Ny - 1; ++i) {
1483  if (thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i + 1].cotbeta) {
1484  ilow = i;
1485  yratio_ = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta) /
1486  (thePixelTemp_[index_id_].enty[i + 1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
1487  break;
1488  }
1489  }
1490  } else {
1491  success_ = false;
1492  }
1493  }
1494 
1495  ihigh = ilow + 1;
1496 
1497  // Use pointers to the three angle pairs used in the interpolation
1498  //
1499  enty0_ = &thePixelTemp_[index_id_].enty[ilow];
1500  enty1_ = &thePixelTemp_[index_id_].enty[ihigh];
1501 
1502  // Interpolate/store all y-related quantities (flip displacements when flip_y_)
1503 
1504  qavg_ = (1.f - yratio_) * enty0_->qavg + yratio_ * enty1_->qavg;
1505  qavg_ *= qcorrect;
1506  symax = (1.f - yratio_) * enty0_->symax + yratio_ * enty1_->symax;
1507  syparmax_ = symax;
1508  sxmax = (1.f - yratio_) * enty0_->sxmax + yratio_ * enty1_->sxmax;
1509  dyone_ = (1.f - yratio_) * enty0_->dyone + yratio_ * enty1_->dyone;
1510  if (flip_y_) {
1511  dyone_ = -dyone_;
1512  }
1513  syone_ = (1.f - yratio_) * enty0_->syone + yratio_ * enty1_->syone;
1514  dytwo_ = (1.f - yratio_) * enty0_->dytwo + yratio_ * enty1_->dytwo;
1515  if (flip_y_) {
1516  dytwo_ = -dytwo_;
1517  }
1518  sytwo_ = (1.f - yratio_) * enty0_->sytwo + yratio_ * enty1_->sytwo;
1519  qmin_ = (1.f - yratio_) * enty0_->qmin + yratio_ * enty1_->qmin;
1520  qmin_ *= qcorrect;
1521  qmin2_ = (1.f - yratio_) * enty0_->qmin2 + yratio_ * enty1_->qmin2;
1522  qmin2_ *= qcorrect;
1523  mpvvav_ = (1.f - yratio_) * enty0_->mpvvav + yratio_ * enty1_->mpvvav;
1524  mpvvav_ *= qcorrect;
1525  sigmavav_ = (1.f - yratio_) * enty0_->sigmavav + yratio_ * enty1_->sigmavav;
1526  kappavav_ = (1.f - yratio_) * enty0_->kappavav + yratio_ * enty1_->kappavav;
1527  mpvvav2_ = (1.f - yratio_) * enty0_->mpvvav2 + yratio_ * enty1_->mpvvav2;
1528  mpvvav2_ *= qcorrect;
1529  sigmavav2_ = (1.f - yratio_) * enty0_->sigmavav2 + yratio_ * enty1_->sigmavav2;
1530  kappavav2_ = (1.f - yratio_) * enty0_->kappavav2 + yratio_ * enty1_->kappavav2;
1531  clsleny_ = fminf(enty0_->clsleny, enty1_->clsleny);
1532  qavg_avg_ = (1.f - yratio_) * enty0_->qavg_avg + yratio_ * enty1_->qavg_avg;
1533  qavg_avg_ *= qcorrect;
1534  for (i = 0; i < 2; ++i) {
1535  for (j = 0; j < 5; ++j) {
1536  // Charge loss switches sides when cot(beta) changes sign
1537 
1538  if (flip_y_) {
1539  yparl_[1 - i][j] = enty0_->ypar[i][j];
1540  yparh_[1 - i][j] = enty1_->ypar[i][j];
1541  } else {
1542  yparl_[i][j] = enty0_->ypar[i][j];
1543  yparh_[i][j] = enty1_->ypar[i][j];
1544  }
1545  if (flip_x_) {
1546  xparly0_[1 - i][j] = enty0_->xpar[i][j];
1547  xparhy0_[1 - i][j] = enty1_->xpar[i][j];
1548  } else {
1549  xparly0_[i][j] = enty0_->xpar[i][j];
1550  xparhy0_[i][j] = enty1_->xpar[i][j];
1551  }
1552  }
1553  }
1554 
1555  for (i = 0; i < 4; ++i) {
1556  yavg_[i] = (1.f - yratio_) * thePixelTemp_[index_id_].enty[ilow].yavg[i] +
1557  yratio_ * thePixelTemp_[index_id_].enty[ihigh].yavg[i];
1558  if (flip_y_) {
1559  yavg_[i] = -yavg_[i];
1560  }
1561  yavg_[i] = (1.f - yratio_) * enty0_->yavg[i] + yratio_ * enty1_->yavg[i];
1562  if (flip_y_) {
1563  yavg_[i] = -yavg_[i];
1564  }
1565  yrms_[i] = (1.f - yratio_) * enty0_->yrms[i] + yratio_ * enty1_->yrms[i];
1566  chi2yavg_[i] = (1.f - yratio_) * enty0_->chi2yavg[i] + yratio_ * enty1_->chi2yavg[i];
1567  chi2ymin_[i] = (1.f - yratio_) * enty0_->chi2ymin[i] + yratio_ * enty1_->chi2ymin[i];
1568  chi2xavg[i] = (1.f - yratio_) * enty0_->chi2xavg[i] + yratio_ * enty1_->chi2xavg[i];
1569  chi2xmin[i] = (1.f - yratio_) * enty0_->chi2xmin[i] + yratio_ * enty1_->chi2xmin[i];
1570  yavgc2m_[i] = (1.f - yratio_) * enty0_->yavgc2m[i] + yratio_ * enty1_->yavgc2m[i];
1571  if (flip_y_) {
1572  yavgc2m_[i] = -yavgc2m_[i];
1573  }
1574  yrmsc2m_[i] = (1.f - yratio_) * enty0_->yrmsc2m[i] + yratio_ * enty1_->yrmsc2m[i];
1575  chi2yavgc2m_[i] = (1.f - yratio_) * enty0_->chi2yavgc2m[i] + yratio_ * enty1_->chi2yavgc2m[i];
1576  // if(flip_y_) {chi2yavgc2m_[i] = -chi2yavgc2m_[i];}
1577  chi2yminc2m_[i] = (1.f - yratio_) * enty0_->chi2yminc2m[i] + yratio_ * enty1_->chi2yminc2m[i];
1578  // xrmsc2m[i]=(1.f - yratio_)*enty0_->xrmsc2m[i] + yratio_*enty1_->xrmsc2m[i];
1579  chi2xavgc2m[i] = (1.f - yratio_) * enty0_->chi2xavgc2m[i] + yratio_ * enty1_->chi2xavgc2m[i];
1580  chi2xminc2m[i] = (1.f - yratio_) * enty0_->chi2xminc2m[i] + yratio_ * enty1_->chi2xminc2m[i];
1581  for (j = 0; j < 6; ++j) {
1582  yflparl_[i][j] = enty0_->yflpar[i][j];
1583  yflparh_[i][j] = enty1_->yflpar[i][j];
1584 
1585  // Since Q_fl is odd under cotbeta, it flips qutomatically, change only even terms
1586 
1587  if (flip_y_ && (j == 0 || j == 2 || j == 4)) {
1588  yflparl_[i][j] = -yflparl_[i][j];
1589  yflparh_[i][j] = -yflparh_[i][j];
1590  }
1591  }
1592  }
1593 
1595 
1596  chi2yavgone_ = (1.f - yratio_) * enty0_->chi2yavgone + yratio_ * enty1_->chi2yavgone;
1597  chi2yminone_ = (1.f - yratio_) * enty0_->chi2yminone + yratio_ * enty1_->chi2yminone;
1598  chi2xavgone = (1.f - yratio_) * enty0_->chi2xavgone + yratio_ * enty1_->chi2xavgone;
1599  chi2xminone = (1.f - yratio_) * enty0_->chi2xminone + yratio_ * enty1_->chi2xminone;
1600 
1601  fracyone_ = (1.f - yratio_) * enty0_->fracyone + yratio_ * enty1_->fracyone;
1602  fracytwo_ = (1.f - yratio_) * enty0_->fracytwo + yratio_ * enty1_->fracytwo;
1603  // If using y-spares
1604  // for(i=0; i<10; ++i) {
1605  // pyspare[i]=(1.f - yratio_)*enty0_->yspare[i] + yratio_*enty1_->yspare[i];
1606  // }
1607 
1608  // Interpolate and build the y-template
1609 
1610  for (i = 0; i < 9; ++i) {
1611  ytemp_[i][0] = 0.f;
1612  ytemp_[i][1] = 0.f;
1613  ytemp_[i][BYM2] = 0.f;
1614  ytemp_[i][BYM1] = 0.f;
1615  for (j = 0; j < TYSIZE; ++j) {
1616  // Flip the basic y-template when the cotbeta is negative
1617 
1618  if (flip_y_) {
1619  ytemp_[8 - i][BYM3 - j] = (1.f - yratio_) * enty0_->ytemp[i][j] + yratio_ * enty1_->ytemp[i][j];
1620  } else {
1621  ytemp_[i][j + 2] = (1.f - yratio_) * enty0_->ytemp[i][j] + yratio_ * enty1_->ytemp[i][j];
1622  }
1623  }
1624  }
1625 
1626  // next, loop over all x-angle entries, first, find relevant y-slices
1627 
1628  iylow = 0;
1629  yxratio = 0.f;
1630 
1631  if (abs_cotb_ >= thePixelTemp_[index_id_].entx[Nyx - 1][0].cotbeta) {
1632  iylow = Nyx - 2;
1633  yxratio = 1.f;
1634 
1635  } else if (abs_cotb_ >= thePixelTemp_[index_id_].entx[0][0].cotbeta) {
1636  for (i = 0; i < Nyx - 1; ++i) {
1637  if (thePixelTemp_[index_id_].entx[i][0].cotbeta <= abs_cotb_ &&
1638  abs_cotb_ < thePixelTemp_[index_id_].entx[i + 1][0].cotbeta) {
1639  iylow = i;
1640  yxratio = (abs_cotb_ - thePixelTemp_[index_id_].entx[i][0].cotbeta) /
1641  (thePixelTemp_[index_id_].entx[i + 1][0].cotbeta - thePixelTemp_[index_id_].entx[i][0].cotbeta);
1642  break;
1643  }
1644  }
1645  }
1646 
1647  iyhigh = iylow + 1;
1648 
1649  ilow = 0;
1650  xxratio = 0.f;
1651 
1652  if (cota >= thePixelTemp_[index_id_].entx[0][Nxx - 1].cotalpha) {
1653  ilow = Nxx - 2;
1654  xxratio = 1.f;
1655  success_ = false;
1656 
1657  } else {
1658  if (cota >= thePixelTemp_[index_id_].entx[0][0].cotalpha) {
1659  for (i = 0; i < Nxx - 1; ++i) {
1660  if (thePixelTemp_[index_id_].entx[0][i].cotalpha <= cota &&
1661  cota < thePixelTemp_[index_id_].entx[0][i + 1].cotalpha) {
1662  ilow = i;
1663  xxratio = (cota - thePixelTemp_[index_id_].entx[0][i].cotalpha) /
1664  (thePixelTemp_[index_id_].entx[0][i + 1].cotalpha - thePixelTemp_[index_id_].entx[0][i].cotalpha);
1665  break;
1666  }
1667  }
1668  } else {
1669  success_ = false;
1670  }
1671  }
1672 
1673  ihigh = ilow + 1;
1674 
1675  // Interpolate/store all x-related quantities
1676 
1677  yxratio_ = yxratio;
1678  xxratio_ = xxratio;
1679 
1680  // sxparmax defines the maximum charge for which the parameters xpar are defined (not rescaled by cotbeta)
1681 
1682  sxparmax_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[imaxx][ilow].sxmax +
1683  xxratio * thePixelTemp_[index_id_].entx[imaxx][ihigh].sxmax;
1684  sxmax_ = sxparmax_;
1685  if (thePixelTemp_[index_id_].entx[imaxx][imidy].sxmax != 0.f) {
1686  sxmax_ = sxmax_ / thePixelTemp_[index_id_].entx[imaxx][imidy].sxmax * sxmax;
1687  }
1688  symax_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[imaxx][ilow].symax +
1689  xxratio * thePixelTemp_[index_id_].entx[imaxx][ihigh].symax;
1690  if (thePixelTemp_[index_id_].entx[imaxx][imidy].symax != 0.f) {
1691  symax_ = symax_ / thePixelTemp_[index_id_].entx[imaxx][imidy].symax * symax;
1692  }
1693  dxone_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[0][ilow].dxone +
1694  xxratio * thePixelTemp_[index_id_].entx[0][ihigh].dxone;
1695  if (flip_x_) {
1696  dxone_ = -dxone_;
1697  }
1698  sxone_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[0][ilow].sxone +
1699  xxratio * thePixelTemp_[index_id_].entx[0][ihigh].sxone;
1700  dxtwo_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[0][ilow].dxtwo +
1701  xxratio * thePixelTemp_[index_id_].entx[0][ihigh].dxtwo;
1702  if (flip_x_) {
1703  dxtwo_ = -dxtwo_;
1704  }
1705  sxtwo_ = (1.f - xxratio) * thePixelTemp_[index_id_].entx[0][ilow].sxtwo +
1706  xxratio * thePixelTemp_[index_id_].entx[0][ihigh].sxtwo;
1707  clslenx_ = fminf(thePixelTemp_[index_id_].entx[0][ilow].clslenx, thePixelTemp_[index_id_].entx[0][ihigh].clslenx);
1708 
1709  for (i = 0; i < 2; ++i) {
1710  for (j = 0; j < 5; ++j) {
1711  // Charge loss switches sides when cot(alpha) changes sign
1712  if (flip_x_) {
1713  xpar0_[1 - i][j] = thePixelTemp_[index_id_].entx[imaxx][imidy].xpar[i][j];
1714  xparl_[1 - i][j] = thePixelTemp_[index_id_].entx[imaxx][ilow].xpar[i][j];
1715  xparh_[1 - i][j] = thePixelTemp_[index_id_].entx[imaxx][ihigh].xpar[i][j];
1716  } else {
1717  xpar0_[i][j] = thePixelTemp_[index_id_].entx[imaxx][imidy].xpar[i][j];
1718  xparl_[i][j] = thePixelTemp_[index_id_].entx[imaxx][ilow].xpar[i][j];
1719  xparh_[i][j] = thePixelTemp_[index_id_].entx[imaxx][ihigh].xpar[i][j];
1720  }
1721  }
1722  }
1723  // Pointers to the currently interpolated point.
1724  entx00_ = &thePixelTemp_[index_id_].entx[iylow][ilow];
1725  entx02_ = &thePixelTemp_[index_id_].entx[iylow][ihigh];
1726  entx20_ = &thePixelTemp_[index_id_].entx[iyhigh][ilow];
1727  entx22_ = &thePixelTemp_[index_id_].entx[iyhigh][ihigh];
1728  entx21_ = &thePixelTemp_[index_id_].entx[iyhigh][imidy];
1729 
1730  // pixmax is the maximum allowed pixel charge (used for truncation)
1731  pixmax_ = (1.f - yxratio) * ((1.f - xxratio) * entx00_->pixmax + xxratio * entx02_->pixmax) +
1732  yxratio * ((1.f - xxratio) * entx20_->pixmax + xxratio * entx22_->pixmax);
1733 
1734  r_qMeas_qTrue_ = (1.f - yxratio) * ((1.f - xxratio) * entx00_->r_qMeas_qTrue + xxratio * entx02_->r_qMeas_qTrue) +
1735  yxratio * ((1.f - xxratio) * entx20_->r_qMeas_qTrue + xxratio * entx22_->r_qMeas_qTrue);
1736 
1737  for (i = 0; i < 4; ++i) {
1738  xavg_[i] = (1.f - yxratio) * ((1.f - xxratio) * entx00_->xavg[i] + xxratio * entx02_->xavg[i]) +
1739  yxratio * ((1.f - xxratio) * entx20_->xavg[i] + xxratio * entx22_->xavg[i]);
1740  if (flip_x_) {
1741  xavg_[i] = -xavg_[i];
1742  }
1743 
1744  xrms_[i] = (1.f - yxratio) * ((1.f - xxratio) * entx00_->xrms[i] + xxratio * entx02_->xrms[i]) +
1745  yxratio * ((1.f - xxratio) * entx20_->xrms[i] + xxratio * entx22_->xrms[i]);
1746 
1747  xavgc2m_[i] = (1.f - yxratio) * ((1.f - xxratio) * entx00_->xavgc2m[i] + xxratio * entx02_->xavgc2m[i]) +
1748  yxratio * ((1.f - xxratio) * entx20_->xavgc2m[i] + xxratio * entx22_->xavgc2m[i]);
1749  if (flip_x_) {
1750  xavgc2m_[i] = -xavgc2m_[i];
1751  }
1752 
1753  xrmsc2m_[i] = (1.f - yxratio) * ((1.f - xxratio) * entx00_->xrmsc2m[i] + xxratio * entx02_->xrmsc2m[i]) +
1754  yxratio * ((1.f - xxratio) * entx20_->xrmsc2m[i] + xxratio * entx22_->xrmsc2m[i]);
1755  //
1756  // Try new interpolation scheme instead
1757  //
1758  //
1759  // chi2xavgc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*entx00_->chi2xavgc2m[i] + xxratio*entx02_->chi2xavgc2m[i])
1760  // +yxratio*((1.f - xxratio)*entx20_->chi2xavgc2m[i] + xxratio*entx22_->chi2xavgc2m[i]);
1761 
1762  // chi2xminc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*entx00_->chi2xminc2m[i] + xxratio*entx02_->chi2xminc2m[i])
1763  // +yxratio*((1.f - xxratio)*entx20_->chi2xminc2m[i] + xxratio*entx22_->chi2xminc2m[i]);
1764  //
1765  // chi2xavg_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].chi2xavg[i] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].chi2xavg[i]);
1766  // if(thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i] != 0.f) {chi2xavg_[i]=chi2xavg_[i]/thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i]*chi2xavg[i];}
1767  //
1768  // chi2xmin_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].chi2xmin[i] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].chi2xmin[i]);
1769  // if(thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i] != 0.f) {chi2xmin_[i]=chi2xmin_[i]/thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i]*chi2xmin[i];}
1770  //
1771  chi2xavg_[i] = ((1.f - xxratio) * entx20_->chi2xavg[i] + xxratio * entx22_->chi2xavg[i]);
1772  if (entx21_->chi2xavg[i] != 0.f) {
1773  chi2xavg_[i] = chi2xavg_[i] / entx21_->chi2xavg[i] * chi2xavg[i];
1774  }
1775 
1776  chi2xmin_[i] = ((1.f - xxratio) * entx20_->chi2xmin[i] + xxratio * entx22_->chi2xmin[i]);
1777  if (entx21_->chi2xmin[i] != 0.f) {
1778  chi2xmin_[i] = chi2xmin_[i] / entx21_->chi2xmin[i] * chi2xmin[i];
1779  }
1780 
1781  chi2xavgc2m_[i] = ((1.f - xxratio) * entx20_->chi2xavgc2m[i] + xxratio * entx22_->chi2xavgc2m[i]);
1782  if (entx21_->chi2xavgc2m[i] != 0.f) {
1783  chi2xavgc2m_[i] = chi2xavgc2m_[i] / entx21_->chi2xavgc2m[i] * chi2xavgc2m[i];
1784  }
1785 
1786  chi2xminc2m_[i] = ((1.f - xxratio) * entx20_->chi2xminc2m[i] + xxratio * entx22_->chi2xminc2m[i]);
1787  if (entx21_->chi2xminc2m[i] != 0.f) {
1788  chi2xminc2m_[i] = chi2xminc2m_[i] / entx21_->chi2xminc2m[i] * chi2xminc2m[i];
1789  }
1790 
1791  for (j = 0; j < 6; ++j) {
1792  xflparll_[i][j] = entx00_->xflpar[i][j];
1793  xflparlh_[i][j] = entx02_->xflpar[i][j];
1794  xflparhl_[i][j] = entx20_->xflpar[i][j];
1795  xflparhh_[i][j] = entx22_->xflpar[i][j];
1796  // Since Q_fl is odd under cotalpha, it flips qutomatically, change only even terms
1797  if (flip_x_ && (j == 0 || j == 2 || j == 4)) {
1798  xflparll_[i][j] = -xflparll_[i][j];
1799  xflparlh_[i][j] = -xflparlh_[i][j];
1800  xflparhl_[i][j] = -xflparhl_[i][j];
1801  xflparhh_[i][j] = -xflparhh_[i][j];
1802  }
1803  }
1804  }
1805 
1806  // Do the spares next
1807 
1808  chi2xavgone_ = ((1.f - xxratio) * entx20_->chi2xavgone + xxratio * entx22_->chi2xavgone);
1809  if (entx21_->chi2xavgone != 0.f) {
1810  chi2xavgone_ = chi2xavgone_ / entx21_->chi2xavgone * chi2xavgone;
1811  }
1812 
1813  chi2xminone_ = ((1.f - xxratio) * entx20_->chi2xminone + xxratio * entx22_->chi2xminone);
1814  if (entx21_->chi2xminone != 0.f) {
1815  chi2xminone_ = chi2xminone_ / entx21_->chi2xminone * chi2xminone;
1816  }
1817 
1818  fracxone_ = (1.f - yxratio) * ((1.f - xxratio) * entx00_->fracxone + xxratio * entx02_->fracxone) +
1819  yxratio * ((1.f - xxratio) * entx20_->fracxone + xxratio * entx22_->fracxone);
1820  fracxtwo_ = (1.f - yxratio) * ((1.f - xxratio) * entx00_->fracxtwo + xxratio * entx02_->fracxtwo) +
1821  yxratio * ((1.f - xxratio) * entx20_->fracxtwo + xxratio * entx22_->fracxtwo);
1822 
1823  // If using x-spares
1824  // for(i=0; i<10; ++i) {
1825  // pxspare[i]=(1.f - yxratio)*((1.f - xxratio)*entx00_->xspare[i] + xxratio*entx02_->xspare[i])
1826  // +yxratio*((1.f - xxratio)*entx20_->xspare[i] + xxratio*entx22_->xspare[i]);
1827  // }
1828 
1829  // Interpolate and build the x-template
1830 
1831  // qxtempcor corrects the total charge to the actual track angles (not actually needed for the template fits, but useful for Guofan)
1832 
1833  cotbeta0 = thePixelTemp_[index_id_].entx[iyhigh][0].cotbeta;
1834  qxtempcor =
1835  std::sqrt((1.f + cotbeta * cotbeta + cotalpha * cotalpha) / (1.f + cotbeta0 * cotbeta0 + cotalpha * cotalpha));
1836 
1837  for (i = 0; i < 9; ++i) {
1838  xtemp_[i][0] = 0.f;
1839  xtemp_[i][1] = 0.f;
1840  xtemp_[i][BXM2] = 0.f;
1841  xtemp_[i][BXM1] = 0.f;
1842  for (j = 0; j < TXSIZE; ++j) {
1843  // Take next largest x-slice for the x-template (it reduces bias in the forward direction after irradiation)
1844  // xtemp_[i][j+2]=(1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].xtemp[i][j] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].xtemp[i][j];
1845  // xtemp_[i][j+2]=(1.f - xxratio)*entx20_->xtemp[i][j] + xxratio*entx22_->xtemp[i][j];
1846  if (flip_x_) {
1847  xtemp_[8 - i][BXM3 - j] =
1848  qxtempcor * ((1.f - xxratio) * entx20_->xtemp[i][j] + xxratio * entx22_->xtemp[i][j]);
1849  } else {
1850  xtemp_[i][j + 2] = qxtempcor * ((1.f - xxratio) * entx20_->xtemp[i][j] + xxratio * entx22_->xtemp[i][j]);
1851  }
1852  }
1853  }
1854 
1855  lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
1856  lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
1857  lorybias_ = thePixelTemp_[index_id_].head.lorybias;
1858  lorxbias_ = thePixelTemp_[index_id_].head.lorxbias;
1859  if (flip_x_) {
1860  lorxwidth_ = -lorxwidth_;
1861  lorxbias_ = -lorxbias_;
1862  }
1863  if (flip_y_) {
1864  lorywidth_ = -lorywidth_;
1865  lorybias_ = -lorybias_;
1866  }
1867  }
1868 
1869  return success_;
1870 } // interpolate
1871 
1872 // ************************************************************************************************************
1878 // ************************************************************************************************************
1879 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta) {
1880  // Interpolate for a new set of track angles
1881 
1882  // Local variables
1883  float locBx = 1.f;
1884  if (cotbeta < 0.f) {
1885  locBx = -1.f;
1886  }
1887  float locBz = locBx;
1888  if (cotalpha < 0.f) {
1889  locBz = -locBx;
1890  }
1891  return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx);
1892 }
1893 
1894 // ************************************************************************************************************
1900 // ************************************************************************************************************
1901 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz) {
1902  // Interpolate for a new set of track angles
1903 
1904  // Local variables
1905  float locBx = 1.f;
1906  return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx);
1907 }
1908 
1909 // *************************************************************************************************************************************
1915 // *************************************************************************************************************************************
1916 #ifdef SI_PIXEL_TEMPLATE_STANDALONE
1917 void SiPixelTemplate::sideload(SiPixelTemplateEntry* entry,
1918  int iDtype,
1919  float locBx,
1920  float locBz,
1921  float lorwdy,
1922  float lorwdx,
1923  float q50,
1924  float fbin[3],
1925  float xsize,
1926  float ysize,
1927  float zsize) {
1928  // Set class variables to the input parameters
1929  entry_sideloaded_ = entry; // will bypass xtemp_[] and ytemp_[] and use those from this Entry.
1930 
1931  enty1_ = entry;
1932  enty0_ = entry;
1933 
1934  entx00_ = entry;
1935  entx02_ = entry;
1936  entx20_ = entry;
1937  entx22_ = entry;
1938  entx21_ = entry;
1939 
1940  dtype_ = iDtype;
1941  lorywidth_ = lorwdy;
1942  lorxwidth_ = lorwdx;
1943  xsize_ = xsize;
1944  ysize_ = ysize;
1945  zsize_ = zsize;
1946  s50_ = q50;
1947  qscale_ = 1.f;
1948  for (int i = 0; i < 3; ++i) {
1949  fbin_[i] = fbin[i];
1950  }
1951 
1952  // Set other class variables
1953 
1954  yratio_ = 0.f;
1955  yxratio_ = 0.f;
1956  xxratio_ = 0.f;
1957 
1958  qavg_ = entry->qavg;
1959  qmin_ = 0.f;
1960  qmin2_ = 0.f;
1961 
1962  pixmax_ = entry->pixmax;
1963  sxmax_ = entry->sxmax;
1964  symax_ = entry->symax;
1965  clsleny_ = entry->clsleny;
1966  clslenx_ = entry->clslenx;
1967 
1968  scaleyavg_ = 1.f;
1969  scalexavg_ = 1.f;
1970  delyavg_ = 0.f;
1971  delysig_ = 0.f;
1972 
1973  dyone_ = entry->dyone;
1974  syone_ = entry->syone;
1975  dytwo_ = entry->dytwo;
1976  sytwo_ = entry->sytwo;
1977 
1978  dxone_ = entry->dxone;
1979  sxone_ = entry->sxone;
1980  dxtwo_ = entry->dxtwo;
1981  sxtwo_ = entry->sxtwo;
1982 
1983  chi2yminone_ = 0.f;
1984  chi2yavgone_ = 0.1;
1985  chi2xminone_ = 0.f;
1986  chi2xavgone_ = 0.1;
1987 
1988  for (int i = 0; i < 4; ++i) {
1989  scalex_[i] = 1.f;
1990  scaley_[i] = 1.f;
1991  offsetx_[i] = 0.f;
1992  offsety_[i] = 0.f;
1993  xrms_[i] = 0.f;
1994  yrms_[i] = 0.f;
1995  xavg_[i] = 0.f;
1996  yavg_[i] = 0.f;
1997 
1998  chi2yavg_[i] = 0.1;
1999  chi2xavg_[i] = 0.1;
2000 
2001  chi2xmin_[i] = 0.f;
2002  chi2ymin_[i] = 0.f;
2003 
2004  for (int j = 0; j < 6; j++) {
2005  yflparl_[i][j] = yflparh_[i][j] = entry->yflpar[i][j];
2006  xflparhh_[i][j] = xflparhl_[i][j] = xflparll_[i][j] = xflparlh_[i][j] = entry->xflpar[i][j];
2007  }
2008  }
2009 
2010  sxparmax_ = entry->sxmax;
2011  syparmax_ = entry->symax;
2012 
2013  // This works only for IP-related tracks
2014 
2015  flip_x_ = false;
2016  flip_y_ = false;
2017  float cotbeta = entry->cotbeta;
2018  switch (dtype_) {
2019  case 0:
2020  if (cotbeta < 0.f) {
2021  flip_y_ = true;
2022  }
2023  break;
2024  case 1:
2025  if (locBz > 0.f) {
2026  flip_y_ = true;
2027  }
2028  break;
2029  case 2:
2030  case 3:
2031  case 4:
2032  case 5:
2033  if (locBx * locBz < 0.f) {
2034  flip_x_ = true;
2035  }
2036  if (locBx < 0.f) {
2037  flip_y_ = true;
2038  }
2039  break;
2040  default:
2041  std::cout << "SiPixelTemplate:illegal subdetector ID = " << dtype_ << std::endl;
2042  }
2043 
2044  // Calculate signed quantities
2045 
2046  // qxtempcor corrects the total charge to the actual track angles (not actually needed for the template fits, but useful for Guofan)
2047  // &&& What to do here?
2048  // qxtempcor=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta0*cotbeta0+cotalpha*cotalpha));
2049  float qxtempcor = 1.f;
2050 
2051  for (int i = 0; i < 9; ++i) {
2052  xtemp_[i][0] = 0.f;
2053  xtemp_[i][1] = 0.f;
2054  xtemp_[i][BXM2] = 0.f;
2055  xtemp_[i][BXM1] = 0.f;
2056  for (int j = 0; j < TXSIZE; ++j) {
2057  if (flip_x_) {
2058  xtemp_[8 - i][BXM3 - j] = qxtempcor * entry_sideloaded_->xtemp[i][j];
2059  } else {
2060  xtemp_[i][j + 2] = qxtempcor * entry_sideloaded_->xtemp[i][j];
2061  }
2062  }
2063  }
2064 
2065  for (int i = 0; i < 9; ++i) {
2066  ytemp_[i][0] = 0.f;
2067  ytemp_[i][1] = 0.f;
2068  ytemp_[i][BYM2] = 0.f;
2069  ytemp_[i][BYM1] = 0.f;
2070  for (int j = 0; j < TYSIZE; ++j) {
2071  // Flip the basic y-template when the cotbeta is negative
2072 
2073  if (flip_y_) {
2074  ytemp_[8 - i][BYM3 - j] = entry_sideloaded_->ytemp[i][j];
2075  } else {
2076  ytemp_[i][j + 2] = entry_sideloaded_->ytemp[i][j];
2077  }
2078  }
2079  }
2080 
2081  // Fitted errors params
2082  for (int i = 0; i < 2; i++) {
2083  for (int j = 0; j < 5; j++) {
2084  if (flip_y_) {
2085  yparl_[1 - i][j] = yparh_[1 - i][j] = entry->ypar[i][j];
2086  } else {
2087  yparl_[i][j] = yparh_[i][j] = entry->ypar[i][j];
2088  }
2089 
2090  if (flip_x_) {
2091  xpar0_[1 - i][j] = xparly0_[1 - i][j] = xparhy0_[1 - i][j] = xparl_[1 - i][j] = xparh_[1 - i][j] =
2092  entry->xpar[i][j];
2093  } else {
2094  xpar0_[i][j] = xparly0_[i][j] = xparhy0_[i][j] = xparl_[i][j] = xparh_[i][j] = entry->xpar[i][j];
2095  }
2096  }
2097  }
2098 
2099  return;
2100 } // sideload
2101 #endif
2102 
2103 // ************************************************************************************************************
2111 // ************************************************************************************************************
2112 void SiPixelTemplate::ysigma2(int fypix, int lypix, float sythr, float ysum[25], float ysig2[25])
2113 
2114 {
2115  // Interpolate using quantities already stored in the private variables
2116 
2117  // Local variables
2118  int i;
2119  float sigi, sigi2, sigi3, sigi4, symax, qscale, s25;
2120 
2121  // Make sure that input is OK
2122 
2123 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2124  if (fypix < 2 || fypix >= BYM2) {
2125  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with fypix = " << fypix << std::endl;
2126  }
2127 #else
2128  assert(fypix > 1 && fypix < BYM2);
2129 #endif
2130 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2131  if (lypix < fypix || lypix >= BYM2) {
2132  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with lypix/fypix = " << lypix << "/"
2133  << fypix << std::endl;
2134  }
2135 #else
2136  assert(lypix >= fypix && lypix < BYM2);
2137 #endif
2138 
2139  // Define the maximum signal to use in the parameterization
2140 
2141  symax = symax_;
2142  s25 = 0.5f * s50_;
2143  if (symax_ > syparmax_) {
2144  symax = syparmax_;
2145  }
2146 
2147  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
2148 
2149  for (i = fypix - 2; i <= lypix + 2; ++i) {
2150  if (i < fypix || i > lypix) {
2151  // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
2152 
2153  ysig2[i] = s50_ * s50_;
2154  } else {
2155  if (ysum[i] < symax) {
2156  sigi = ysum[i];
2157  qscale = 1.f;
2158  if (sigi < s25)
2159  sigi = s25;
2160  } else {
2161  sigi = symax;
2162  qscale = ysum[i] / symax;
2163  }
2164  sigi2 = sigi * sigi;
2165  sigi3 = sigi2 * sigi;
2166  sigi4 = sigi3 * sigi;
2167  if (i <= BHY) {
2168  ysig2[i] = (1.f - yratio_) * (yparl_[0][0] + yparl_[0][1] * sigi + yparl_[0][2] * sigi2 + yparl_[0][3] * sigi3 +
2169  yparl_[0][4] * sigi4) +
2170  yratio_ * (yparh_[0][0] + yparh_[0][1] * sigi + yparh_[0][2] * sigi2 + yparh_[0][3] * sigi3 +
2171  yparh_[0][4] * sigi4);
2172  } else {
2173  ysig2[i] = (1.f - yratio_) * (yparl_[1][0] + yparl_[1][1] * sigi + yparl_[1][2] * sigi2 + yparl_[1][3] * sigi3 +
2174  yparl_[1][4] * sigi4) +
2175  yratio_ * (yparh_[1][0] + yparh_[1][1] * sigi + yparh_[1][2] * sigi2 + yparh_[1][3] * sigi3 +
2176  yparh_[1][4] * sigi4);
2177  }
2178  ysig2[i] *= qscale;
2179  if (ysum[i] > sythr) {
2180  ysig2[i] = 1.e8;
2181  }
2182  if (ysig2[i] <= 0.f) {
2183  LOGERROR("SiPixelTemplate") << "neg y-error-squared, id = " << id_current_ << ", index = " << index_id_
2184  << ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_
2185  << ", sigi = " << sigi << ENDL;
2186  }
2187  }
2188  }
2189 
2190  return;
2191 
2192 } // End ysigma2
2193 
2194 // ************************************************************************************************************
2200 // ************************************************************************************************************
2201 void SiPixelTemplate::ysigma2(float qpixel, int index, float& ysig2)
2202 
2203 {
2204  // Interpolate using quantities already stored in the private variables
2205 
2206  // Local variables
2207  float sigi, sigi2, sigi3, sigi4, symax, qscale, err2, s25;
2208 
2209  // Make sure that input is OK
2210 
2211 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2212  if (index < 2 || index >= BYM2) {
2213  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with index = " << index << std::endl;
2214  }
2215 #else
2216  assert(index > 1 && index < BYM2);
2217 #endif
2218 
2219  // Define the maximum signal to use in the parameterization
2220 
2221  symax = symax_;
2222  s25 = 0.5f * s50_;
2223  if (symax_ > syparmax_) {
2224  symax = syparmax_;
2225  }
2226 
2227  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
2228 
2229  if (qpixel < symax) {
2230  sigi = qpixel;
2231  qscale = 1.f;
2232  if (sigi < s25)
2233  sigi = s25;
2234  } else {
2235  sigi = symax;
2236  qscale = qpixel / symax;
2237  }
2238  sigi2 = sigi * sigi;
2239  sigi3 = sigi2 * sigi;
2240  sigi4 = sigi3 * sigi;
2241  if (index <= BHY) {
2242  err2 =
2243  (1.f - yratio_) *
2244  (yparl_[0][0] + yparl_[0][1] * sigi + yparl_[0][2] * sigi2 + yparl_[0][3] * sigi3 + yparl_[0][4] * sigi4) +
2245  yratio_ *
2246  (yparh_[0][0] + yparh_[0][1] * sigi + yparh_[0][2] * sigi2 + yparh_[0][3] * sigi3 + yparh_[0][4] * sigi4);
2247  } else {
2248  err2 =
2249  (1.f - yratio_) *
2250  (yparl_[1][0] + yparl_[1][1] * sigi + yparl_[1][2] * sigi2 + yparl_[1][3] * sigi3 + yparl_[1][4] * sigi4) +
2251  yratio_ *
2252  (yparh_[1][0] + yparh_[1][1] * sigi + yparh_[1][2] * sigi2 + yparh_[1][3] * sigi3 + yparh_[1][4] * sigi4);
2253  }
2254  ysig2 = qscale * err2;
2255  if (ysig2 <= 0.f) {
2256  LOGERROR("SiPixelTemplate") << "neg y-error-squared, id = " << id_current_ << ", index = " << index_id_
2257  << ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_
2258  << ", sigi = " << sigi << ENDL;
2259  }
2260 
2261  return;
2262 
2263 } // End ysigma2
2264 
2265 // ************************************************************************************************************
2273 // ************************************************************************************************************
2274 void SiPixelTemplate::xsigma2(int fxpix, int lxpix, float sxthr, float xsum[11], float xsig2[11])
2275 
2276 {
2277  // Interpolate using quantities already stored in the private variables
2278 
2279  // Local variables
2280  int i;
2281  float sigi, sigi2, sigi3, sigi4, yint, sxmax, x0, qscale, s25;
2282 
2283  // Make sure that input is OK
2284 
2285 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2286  if (fxpix < 2 || fxpix >= BXM2) {
2287  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xsigma2 called with fxpix = " << fxpix << std::endl;
2288  }
2289 #else
2290  assert(fxpix > 1 && fxpix < BXM2);
2291 #endif
2292 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2293  if (lxpix < fxpix || lxpix >= BXM2) {
2294  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xsigma2 called with lxpix/fxpix = " << lxpix << "/"
2295  << fxpix << std::endl;
2296  }
2297 #else
2298  assert(lxpix >= fxpix && lxpix < BXM2);
2299 #endif
2300 
2301  // Define the maximum signal to use in the parameterization
2302 
2303  sxmax = sxmax_;
2304  s25 = 0.5f * s50_;
2305  if (sxmax_ > sxparmax_) {
2306  sxmax = sxparmax_;
2307  }
2308 
2309  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
2310 
2311  for (i = fxpix - 2; i <= lxpix + 2; ++i) {
2312  if (i < fxpix || i > lxpix) {
2313  // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
2314 
2315  xsig2[i] = s50_ * s50_;
2316  } else {
2317  if (xsum[i] < sxmax) {
2318  sigi = xsum[i];
2319  qscale = 1.f;
2320  if (sigi < s25)
2321  sigi = s25;
2322  } else {
2323  sigi = sxmax;
2324  qscale = xsum[i] / sxmax;
2325  }
2326  sigi2 = sigi * sigi;
2327  sigi3 = sigi2 * sigi;
2328  sigi4 = sigi3 * sigi;
2329 
2330  // First, do the cotbeta interpolation
2331 
2332  if (i <= BHX) {
2333  yint = (1.f - yratio_) * (xparly0_[0][0] + xparly0_[0][1] * sigi + xparly0_[0][2] * sigi2 +
2334  xparly0_[0][3] * sigi3 + xparly0_[0][4] * sigi4) +
2335  yratio_ * (xparhy0_[0][0] + xparhy0_[0][1] * sigi + xparhy0_[0][2] * sigi2 + xparhy0_[0][3] * sigi3 +
2336  xparhy0_[0][4] * sigi4);
2337  } else {
2338  yint = (1.f - yratio_) * (xparly0_[1][0] + xparly0_[1][1] * sigi + xparly0_[1][2] * sigi2 +
2339  xparly0_[1][3] * sigi3 + xparly0_[1][4] * sigi4) +
2340  yratio_ * (xparhy0_[1][0] + xparhy0_[1][1] * sigi + xparhy0_[1][2] * sigi2 + xparhy0_[1][3] * sigi3 +
2341  xparhy0_[1][4] * sigi4);
2342  }
2343 
2344  // Next, do the cotalpha interpolation
2345 
2346  if (i <= BHX) {
2347  xsig2[i] = (1.f - xxratio_) * (xparl_[0][0] + xparl_[0][1] * sigi + xparl_[0][2] * sigi2 +
2348  xparl_[0][3] * sigi3 + xparl_[0][4] * sigi4) +
2349  xxratio_ * (xparh_[0][0] + xparh_[0][1] * sigi + xparh_[0][2] * sigi2 + xparh_[0][3] * sigi3 +
2350  xparh_[0][4] * sigi4);
2351  } else {
2352  xsig2[i] = (1.f - xxratio_) * (xparl_[1][0] + xparl_[1][1] * sigi + xparl_[1][2] * sigi2 +
2353  xparl_[1][3] * sigi3 + xparl_[1][4] * sigi4) +
2354  xxratio_ * (xparh_[1][0] + xparh_[1][1] * sigi + xparh_[1][2] * sigi2 + xparh_[1][3] * sigi3 +
2355  xparh_[1][4] * sigi4);
2356  }
2357 
2358  // Finally, get the mid-point value of the cotalpha function
2359 
2360  if (i <= BHX) {
2361  x0 = xpar0_[0][0] + xpar0_[0][1] * sigi + xpar0_[0][2] * sigi2 + xpar0_[0][3] * sigi3 + xpar0_[0][4] * sigi4;
2362  } else {
2363  x0 = xpar0_[1][0] + xpar0_[1][1] * sigi + xpar0_[1][2] * sigi2 + xpar0_[1][3] * sigi3 + xpar0_[1][4] * sigi4;
2364  }
2365 
2366  // Finally, rescale the yint value for cotalpha variation
2367 
2368  if (x0 != 0.f) {
2369  xsig2[i] = xsig2[i] / x0 * yint;
2370  }
2371  xsig2[i] *= qscale;
2372  if (xsum[i] > sxthr) {
2373  xsig2[i] = 1.e8;
2374  }
2375  if (xsig2[i] <= 0.f) {
2376  LOGERROR("SiPixelTemplate") << "neg x-error-squared, id = " << id_current_ << ", index = " << index_id_
2377  << ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_
2378  << ", sigi = " << sigi << ENDL;
2379  }
2380  }
2381  }
2382 
2383  return;
2384 
2385 } // End xsigma2
2386 
2387 // ************************************************************************************************************
2391 // ************************************************************************************************************
2392 float SiPixelTemplate::yflcorr(int binq, float qfly)
2393 
2394 {
2395  // Interpolate using quantities already stored in the private variables
2396 
2397  // Local variables
2398  float qfl, qfl2, qfl3, qfl4, qfl5, dy;
2399 
2400  // Make sure that input is OK
2401 
2402 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2403  if (binq < 0 || binq > 3) {
2404  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yflcorr called with binq = " << binq << std::endl;
2405  }
2406 #else
2407  assert(binq >= 0 && binq < 4);
2408 #endif
2409 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2410  if (fabs((double)qfly) > 1.) {
2411  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yflcorr called with qfly = " << qfly << std::endl;
2412  }
2413 #else
2414  assert(fabs((double)qfly) <= 1.);
2415 #endif
2416 
2417  // Define the maximum signal to allow before de-weighting a pixel
2418 
2419  qfl = qfly;
2420 
2421  if (qfl < -0.9f) {
2422  qfl = -0.9f;
2423  }
2424  if (qfl > 0.9f) {
2425  qfl = 0.9f;
2426  }
2427 
2428  // Interpolate between the two polynomials
2429 
2430  qfl2 = qfl * qfl;
2431  qfl3 = qfl2 * qfl;
2432  qfl4 = qfl3 * qfl;
2433  qfl5 = qfl4 * qfl;
2434  dy = (1.f - yratio_) * (yflparl_[binq][0] + yflparl_[binq][1] * qfl + yflparl_[binq][2] * qfl2 +
2435  yflparl_[binq][3] * qfl3 + yflparl_[binq][4] * qfl4 + yflparl_[binq][5] * qfl5) +
2436  yratio_ * (yflparh_[binq][0] + yflparh_[binq][1] * qfl + yflparh_[binq][2] * qfl2 + yflparh_[binq][3] * qfl3 +
2437  yflparh_[binq][4] * qfl4 + yflparh_[binq][5] * qfl5);
2438 
2439  return dy;
2440 
2441 } // End yflcorr
2442 
2443 // ************************************************************************************************************
2447 // ************************************************************************************************************
2448 float SiPixelTemplate::xflcorr(int binq, float qflx)
2449 
2450 {
2451  // Interpolate using quantities already stored in the private variables
2452 
2453  // Local variables
2454  float qfl, qfl2, qfl3, qfl4, qfl5, dx;
2455 
2456  // Make sure that input is OK
2457 
2458 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2459  if (binq < 0 || binq > 3) {
2460  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xflcorr called with binq = " << binq << std::endl;
2461  }
2462 #else
2463  assert(binq >= 0 && binq < 4);
2464 #endif
2465 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2466  if (fabs((double)qflx) > 1.) {
2467  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xflcorr called with qflx = " << qflx << std::endl;
2468  }
2469 #else
2470  assert(fabs((double)qflx) <= 1.);
2471 #endif
2472 
2473  // Define the maximum signal to allow before de-weighting a pixel
2474 
2475  qfl = qflx;
2476 
2477  if (qfl < -0.9f) {
2478  qfl = -0.9f;
2479  }
2480  if (qfl > 0.9f) {
2481  qfl = 0.9f;
2482  }
2483 
2484  // Interpolate between the two polynomials
2485 
2486  qfl2 = qfl * qfl;
2487  qfl3 = qfl2 * qfl;
2488  qfl4 = qfl3 * qfl;
2489  qfl5 = qfl4 * qfl;
2490  dx = (1.f - yxratio_) *
2491  ((1.f - xxratio_) * (xflparll_[binq][0] + xflparll_[binq][1] * qfl + xflparll_[binq][2] * qfl2 +
2492  xflparll_[binq][3] * qfl3 + xflparll_[binq][4] * qfl4 + xflparll_[binq][5] * qfl5) +
2493  xxratio_ * (xflparlh_[binq][0] + xflparlh_[binq][1] * qfl + xflparlh_[binq][2] * qfl2 +
2494  xflparlh_[binq][3] * qfl3 + xflparlh_[binq][4] * qfl4 + xflparlh_[binq][5] * qfl5)) +
2495  yxratio_ *
2496  ((1.f - xxratio_) * (xflparhl_[binq][0] + xflparhl_[binq][1] * qfl + xflparhl_[binq][2] * qfl2 +
2497  xflparhl_[binq][3] * qfl3 + xflparhl_[binq][4] * qfl4 + xflparhl_[binq][5] * qfl5) +
2498  xxratio_ * (xflparhh_[binq][0] + xflparhh_[binq][1] * qfl + xflparhh_[binq][2] * qfl2 +
2499  xflparhh_[binq][3] * qfl3 + xflparhh_[binq][4] * qfl4 + xflparhh_[binq][5] * qfl5));
2500 
2501  return dx;
2502 
2503 } // End xflcorr
2504 
2505 // ************************************************************************************************************
2510 // ************************************************************************************************************
2511 void SiPixelTemplate::ytemp(int fybin, int lybin, float ytemplate[41][BYSIZE])
2512 
2513 {
2514  // Retrieve already interpolated quantities
2515 
2516  // Local variables
2517  int i, j;
2518 
2519  // Verify that input parameters are in valid range
2520 
2521 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2522  if (fybin < 0 || fybin > 40) {
2523  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp called with fybin = " << fybin << std::endl;
2524  }
2525 #else
2526  assert(fybin >= 0 && fybin < 41);
2527 #endif
2528 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2529  if (lybin < 0 || lybin > 40) {
2530  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp called with lybin = " << lybin << std::endl;
2531  }
2532 #else
2533  assert(lybin >= 0 && lybin < 41);
2534 #endif
2535 
2536  // Build the y-template, the central 25 bins are here in all cases
2537 
2538  for (i = 0; i < 9; ++i) {
2539  for (j = 0; j < BYSIZE; ++j) {
2540  ytemplate[i + 16][j] = ytemp_[i][j];
2541  }
2542  }
2543  for (i = 0; i < 8; ++i) {
2544  ytemplate[i + 8][BYM1] = 0.f;
2545  for (j = 0; j < BYM1; ++j) {
2546  ytemplate[i + 8][j] = ytemp_[i][j + 1];
2547  }
2548  }
2549  for (i = 1; i < 9; ++i) {
2550  ytemplate[i + 24][0] = 0.f;
2551  for (j = 0; j < BYM1; ++j) {
2552  ytemplate[i + 24][j + 1] = ytemp_[i][j];
2553  }
2554  }
2555 
2556  // Add more bins if needed
2557 
2558  if (fybin < 8) {
2559  for (i = 0; i < 8; ++i) {
2560  ytemplate[i][BYM2] = 0.f;
2561  ytemplate[i][BYM1] = 0.f;
2562  for (j = 0; j < BYM2; ++j) {
2563  ytemplate[i][j] = ytemp_[i][j + 2];
2564  }
2565  }
2566  }
2567  if (lybin > 32) {
2568  for (i = 1; i < 9; ++i) {
2569  ytemplate[i + 32][0] = 0.f;
2570  ytemplate[i + 32][1] = 0.f;
2571  for (j = 0; j < BYM2; ++j) {
2572  ytemplate[i + 32][j + 2] = ytemp_[i][j];
2573  }
2574  }
2575  }
2576 
2577  return;
2578 
2579 } // End ytemp
2580 
2581 // ************************************************************************************************************
2586 // ************************************************************************************************************
2587 void SiPixelTemplate::xtemp(int fxbin, int lxbin, float xtemplate[41][BXSIZE]) {
2588  // Retrieve already interpolated quantities
2589 
2590  // Local variables
2591  int i, j;
2592 
2593  // Verify that input parameters are in valid range
2594 
2595 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2596  if (fxbin < 0 || fxbin > 40) {
2597  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp called with fxbin = " << fxbin << std::endl;
2598  }
2599 #else
2600  assert(fxbin >= 0 && fxbin < 41);
2601 #endif
2602 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2603  if (lxbin < 0 || lxbin > 40) {
2604  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp called with lxbin = " << lxbin << std::endl;
2605  }
2606 #else
2607  assert(lxbin >= 0 && lxbin < 41);
2608 #endif
2609 
2610  // Build the x-template, the central 25 bins are here in all cases
2611 
2612  for (i = 0; i < 9; ++i) {
2613  for (j = 0; j < BXSIZE; ++j) {
2614  xtemplate[i + 16][j] = xtemp_[i][j];
2615  }
2616  }
2617  for (i = 0; i < 8; ++i) {
2618  xtemplate[i + 8][BXM1] = 0.f;
2619  for (j = 0; j < BXM1; ++j) {
2620  xtemplate[i + 8][j] = xtemp_[i][j + 1];
2621  }
2622  }
2623  for (i = 1; i < 9; ++i) {
2624  xtemplate[i + 24][0] = 0.f;
2625  for (j = 0; j < BXM1; ++j) {
2626  xtemplate[i + 24][j + 1] = xtemp_[i][j];
2627  }
2628  }
2629  // Add more bins if needed
2630  if (fxbin < 8) {
2631  for (i = 0; i < 8; ++i) {
2632  xtemplate[i][BXM2] = 0.f;
2633  xtemplate[i][BXM1] = 0.f;
2634  for (j = 0; j < BXM2; ++j) {
2635  xtemplate[i][j] = xtemp_[i][j + 2];
2636  }
2637  }
2638  }
2639  if (lxbin > 32) {
2640  for (i = 1; i < 9; ++i) {
2641  xtemplate[i + 32][0] = 0.f;
2642  xtemplate[i + 32][1] = 0.f;
2643  for (j = 0; j < BXM2; ++j) {
2644  xtemplate[i + 32][j + 2] = xtemp_[i][j];
2645  }
2646  }
2647  }
2648 
2649  return;
2650 
2651 } // End xtemp
2652 
2653 // ************************************************************************************************************
2655 // ************************************************************************************************************
2656 int SiPixelTemplate::cytemp()
2657 
2658 {
2659  // Retrieve already interpolated quantities
2660 
2661  // Local variables
2662  int j;
2663 
2664  // Analyze only pixels along the central entry
2665  // First, find the maximum signal and then work out to the edges
2666 
2667  float sigmax = 0.f;
2668  float qedge = 2. * s50_;
2669  int jmax = -1;
2670 
2671  for (j = 0; j < BYSIZE; ++j) {
2672  if (ytemp_[4][j] > sigmax) {
2673  sigmax = ytemp_[4][j];
2674  jmax = j;
2675  }
2676  }
2677  if (sigmax < qedge) {
2678  qedge = s50_;
2679  }
2680  if (sigmax < qedge || jmax < 1 || jmax > BYM2) {
2681  return -1;
2682  }
2683 
2684  // Now search forward and backward
2685 
2686  int jend = jmax;
2687 
2688  for (j = jmax + 1; j < BYM1; ++j) {
2689  if (ytemp_[4][j] < qedge)
2690  break;
2691  jend = j;
2692  }
2693 
2694  int jbeg = jmax;
2695 
2696  for (j = jmax - 1; j > 0; --j) {
2697  if (ytemp_[4][j] < qedge)
2698  break;
2699  jbeg = j;
2700  }
2701 
2702  return (jbeg + jend) / 2;
2703 
2704 } // End cytemp
2705 
2706 // ************************************************************************************************************
2708 // ************************************************************************************************************
2709 int SiPixelTemplate::cxtemp()
2710 
2711 {
2712  // Retrieve already interpolated quantities
2713 
2714  // Local variables
2715  int j;
2716 
2717  // Analyze only pixels along the central entry
2718  // First, find the maximum signal and then work out to the edges
2719 
2720  float sigmax = 0.f;
2721  float qedge = 2. * s50_;
2722  int jmax = -1;
2723 
2724  for (j = 0; j < BXSIZE; ++j) {
2725  if (xtemp_[4][j] > sigmax) {
2726  sigmax = xtemp_[4][j];
2727  jmax = j;
2728  }
2729  }
2730  if (sigmax < qedge) {
2731  qedge = s50_;
2732  }
2733  if (sigmax < qedge || jmax < 1 || jmax > BXM2) {
2734  return -1;
2735  }
2736 
2737  // Now search forward and backward
2738 
2739  int jend = jmax;
2740 
2741  for (j = jmax + 1; j < BXM1; ++j) {
2742  if (xtemp_[4][j] < qedge)
2743  break;
2744  jend = j;
2745  }
2746 
2747  int jbeg = jmax;
2748 
2749  for (j = jmax - 1; j > 0; --j) {
2750  if (xtemp_[4][j] < qedge)
2751  break;
2752  jbeg = j;
2753  }
2754 
2755  return (jbeg + jend) / 2;
2756 
2757 } // End cxtemp
2758 
2759 // ************************************************************************************************************
2763 // ************************************************************************************************************
2764 void SiPixelTemplate::ytemp3d_int(int nypix, int& nybins)
2765 
2766 {
2767  // Retrieve already interpolated quantities
2768 
2769  // Local variables
2770  int i, j, k;
2771  int ioff0, ioffp, ioffm;
2772 
2773  // Verify that input parameters are in valid range
2774 
2775 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2776  if (nypix < 1 || nypix >= BYM3) {
2777  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp3d called with nypix = " << nypix << std::endl;
2778  }
2779 #else
2780  assert(nypix > 0 && nypix < BYM3);
2781 #endif
2782 
2783  // Calculate the size of the shift in pixels needed to span the entire cluster
2784 
2785  float diff = fabsf(nypix - clsleny_) / 2. + 1.f;
2786  int nshift = (int)diff;
2787  if ((diff - nshift) > 0.5f) {
2788  ++nshift;
2789  }
2790 
2791  // Calculate the number of bins needed to specify each hit range
2792 
2793  nybins_ = 9 + 16 * nshift;
2794 
2795  // Create a 2-d working template with the correct size
2796 
2797  temp2dy_.resize(boost::extents[nybins_][BYSIZE]);
2798 
2799  // The 9 central bins are copied from the interpolated private store
2800 
2801  ioff0 = 8 * nshift;
2802 
2803  for (i = 0; i < 9; ++i) {
2804  for (j = 0; j < BYSIZE; ++j) {
2805  temp2dy_[i + ioff0][j] = ytemp_[i][j];
2806  }
2807  }
2808 
2809  // Add the +- shifted templates
2810 
2811  for (k = 1; k <= nshift; ++k) {
2812  ioffm = ioff0 - k * 8;
2813  for (i = 0; i < 8; ++i) {
2814  for (j = 0; j < k; ++j) {
2815  temp2dy_[i + ioffm][BYM1 - j] = 0.f;
2816  }
2817  for (j = 0; j < BYSIZE - k; ++j) {
2818  temp2dy_[i + ioffm][j] = ytemp_[i][j + k];
2819  }
2820  }
2821  ioffp = ioff0 + k * 8;
2822  for (i = 1; i < 9; ++i) {
2823  for (j = 0; j < k; ++j) {
2824  temp2dy_[i + ioffp][j] = 0.f;
2825  }
2826  for (j = 0; j < BYSIZE - k; ++j) {
2827  temp2dy_[i + ioffp][j + k] = ytemp_[i][j];
2828  }
2829  }
2830  }
2831 
2832  nybins = nybins_;
2833  return;
2834 
2835 } // End ytemp3d_int
2836 
2837 // ************************************************************************************************************
2841 // ************************************************************************************************************
2842 void SiPixelTemplate::ytemp3d(int i, int j, std::vector<float>& ytemplate)
2843 
2844 {
2845  // Sum two 2-d templates to make the 3-d template
2846  if (i >= 0 && i < nybins_ && j <= i) {
2847  for (int k = 0; k < BYSIZE; ++k) {
2848  ytemplate[k] = temp2dy_[i][k] + temp2dy_[j][k];
2849  }
2850  } else {
2851  for (int k = 0; k < BYSIZE; ++k) {
2852  ytemplate[k] = 0.;
2853  }
2854  }
2855 
2856  return;
2857 
2858 } // End ytemp3d
2859 
2860 // ************************************************************************************************************
2864 // ************************************************************************************************************
2865 void SiPixelTemplate::xtemp3d_int(int nxpix, int& nxbins)
2866 
2867 {
2868  // Retrieve already interpolated quantities
2869 
2870  // Local variables
2871  int i, j, k;
2872  int ioff0, ioffp, ioffm;
2873 
2874  // Verify that input parameters are in valid range
2875 
2876 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2877  if (nxpix < 1 || nxpix >= BXM3) {
2878  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp3d called with nxpix = " << nxpix << std::endl;
2879  }
2880 #else
2881  assert(nxpix > 0 && nxpix < BXM3);
2882 #endif
2883 
2884  // Calculate the size of the shift in pixels needed to span the entire cluster
2885 
2886  float diff = fabsf(nxpix - clslenx_) / 2.f + 1.f;
2887  int nshift = (int)diff;
2888  if ((diff - nshift) > 0.5f) {
2889  ++nshift;
2890  }
2891 
2892  // Calculate the number of bins needed to specify each hit range
2893 
2894  nxbins_ = 9 + 16 * nshift;
2895 
2896  // Create a 2-d working template with the correct size
2897 
2898  temp2dx_.resize(boost::extents[nxbins_][BXSIZE]);
2899 
2900  // The 9 central bins are copied from the interpolated private store
2901 
2902  ioff0 = 8 * nshift;
2903 
2904  for (i = 0; i < 9; ++i) {
2905  for (j = 0; j < BXSIZE; ++j) {
2906  temp2dx_[i + ioff0][j] = xtemp_[i][j];
2907  }
2908  }
2909 
2910  // Add the +- shifted templates
2911 
2912  for (k = 1; k <= nshift; ++k) {
2913  ioffm = ioff0 - k * 8;
2914  for (i = 0; i < 8; ++i) {
2915  for (j = 0; j < k; ++j) {
2916  temp2dx_[i + ioffm][BXM1 - j] = 0.f;
2917  }
2918  for (j = 0; j < BXSIZE - k; ++j) {
2919  temp2dx_[i + ioffm][j] = xtemp_[i][j + k];
2920  }
2921  }
2922  ioffp = ioff0 + k * 8;
2923  for (i = 1; i < 9; ++i) {
2924  for (j = 0; j < k; ++j) {
2925  temp2dx_[i + ioffp][j] = 0.f;
2926  }
2927  for (j = 0; j < BXSIZE - k; ++j) {
2928  temp2dx_[i + ioffp][j + k] = xtemp_[i][j];
2929  }
2930  }
2931  }
2932 
2933  nxbins = nxbins_;
2934 
2935  return;
2936 
2937 } // End xtemp3d_int
2938 
2939 // ************************************************************************************************************
2943 // ************************************************************************************************************
2944 void SiPixelTemplate::xtemp3d(int i, int j, std::vector<float>& xtemplate)
2945 
2946 {
2947  // Sum two 2-d templates to make the 3-d template
2948  if (i >= 0 && i < nxbins_ && j <= i) {
2949  for (int k = 0; k < BXSIZE; ++k) {
2950  xtemplate[k] = temp2dx_[i][k] + temp2dx_[j][k];
2951  }
2952  } else {
2953  for (int k = 0; k < BXSIZE; ++k) {
2954  xtemplate[k] = 0.;
2955  }
2956  }
2957 
2958  return;
2959 
2960 } // End xtemp3d
2961 
2962 // ************************************************************************************************************
2990 // ************************************************************************************************************
2991 int SiPixelTemplate::qbin(int id,
2992  float cotalpha,
2993  float cotbeta,
2994  float locBz,
2995  float locBx,
2996  float qclus,
2997  float& pixmx,
2998  float& sigmay,
2999  float& deltay,
3000  float& sigmax,
3001  float& deltax,
3002  float& sy1,
3003  float& dy1,
3004  float& sy2,
3005  float& dy2,
3006  float& sx1,
3007  float& dx1,
3008  float& sx2,
3009  float& dx2) // float& lorywidth, float& lorxwidth)
3010 {
3011  // Interpolate for a new set of track angles
3012 
3013  // &&& New approach: use cached pointers.
3014 
3015  // Find the index corresponding to id
3016 
3017  int index = -1;
3018  for (int i = 0; i < (int)thePixelTemp_.size(); ++i) {
3019  if (id == thePixelTemp_[i].head.ID) {
3020  index = i;
3021  break;
3022  }
3023  }
3024 
3025 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3026  if (index < 0 || index >= (int)thePixelTemp_.size()) {
3027  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::qbin can't find needed template ID = " << id << std::endl;
3028  }
3029 #else
3030  assert(index >= 0 && index < (int)thePixelTemp_.size());
3031 #endif
3032 
3033  //
3034 
3035  auto const& templ = thePixelTemp_[index];
3036 
3037  // Interpolate the absolute value of cot(beta)
3038 
3039  auto acotb = std::abs(cotbeta);
3040 
3041  // qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
3042 
3043  auto cotalpha0 = thePixelTemp_[index].enty[0].cotalpha;
3044  auto qcorrect =
3045  std::sqrt((1.f + cotbeta * cotbeta + cotalpha * cotalpha) / (1.f + cotbeta * cotbeta + cotalpha0 * cotalpha0));
3046 
3047  // for some cosmics, the ususal gymnastics are incorrect
3048 
3049  float cota = cotalpha;
3050  bool flip_x = false;
3051  //y flipping already taken care of by interpolate
3052  switch (dtype_) {
3053  case 0:
3054  case 1:
3055  case 2:
3056  case 3:
3057  case 4:
3058  case 5:
3059  if (locBx * locBz < 0.f) {
3060  cota = -cotalpha;
3061  flip_x = true;
3062  }
3063  break;
3064  default:
3065 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3066  throw cms::Exception("DataCorrupt")
3067  << "SiPixelTemplate::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
3068 #else
3069  std::cout << "SiPixelTemplate::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
3070 #endif
3071  }
3072 
3073  // Copy the charge scaling factor to the private variable
3074 
3075  auto qscale = thePixelTemp_[index].head.qscale;
3076 
3077  /*
3078  lorywidth = thePixelTemp_[index].head.lorywidth;
3079  if(locBz > 0.f) {lorywidth = -lorywidth;}
3080  lorxwidth = thePixelTemp_[index].head.lorxwidth;
3081  */
3082 
3083  auto Ny = thePixelTemp_[index].head.NTy;
3084  auto Nyx = thePixelTemp_[index].head.NTyx;
3085  auto Nxx = thePixelTemp_[index].head.NTxx;
3086 
3087 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3088  if (Ny < 2 || Nyx < 1 || Nxx < 2) {
3089  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny
3090  << "/" << Nyx << "/" << Nxx << std::endl;
3091  }
3092 #else
3093  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
3094 #endif
3095 
3096  // Interpolate/store all y-related quantities (flip displacements when flip_y)
3097 
3098  dy1 = (1.f - yratio_) * enty0_->dyone + yratio_ * enty1_->dyone;
3099  if (flip_y_) {
3100  dy1 = -dy1;
3101  }
3102  sy1 = (1.f - yratio_) * enty0_->syone + yratio_ * enty1_->syone;
3103  dy2 = (1.f - yratio_) * enty0_->dytwo + yratio_ * enty1_->dytwo;
3104  if (flip_y_) {
3105  dy2 = -dy2;
3106  }
3107  sy2 = (1.f - yratio_) * enty0_->sytwo + yratio_ * enty1_->sytwo;
3108 
3109  auto qavg = (1.f - yratio_) * enty0_->qavg + yratio_ * enty1_->qavg;
3110  qavg *= qcorrect;
3111  auto qmin = (1.f - yratio_) * enty0_->qmin + yratio_ * enty1_->qmin;
3112  qmin *= qcorrect;
3113  auto qmin2 = (1.f - yratio_) * enty0_->qmin2 + yratio_ * enty1_->qmin2;
3114  qmin2 *= qcorrect;
3115 
3116 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3117  if (qavg <= 0.f || qmin <= 0.f) {
3118  throw cms::Exception("DataCorrupt")
3119  << "SiPixelTemplate::qbin, qavg or qmin <= 0,"
3120  << " Probably someone called the generic pixel reconstruction with an illegal trajectory state" << std::endl;
3121  }
3122 #else
3123  assert(qavg > 0.f && qmin > 0.f);
3124 #endif
3125 
3126  // Scale the input charge to account for differences between pixelav and CMSSW simulation or data
3127 
3128  auto qtotal = qscale * qclus;
3129 
3130  // uncertainty and final corrections depend upon total charge bin
3131  auto fq = qtotal / qavg;
3132  int binq;
3133 
3134  if (fq > fbin_[0]) {
3135  binq = 0;
3136  //std::cout<<" fq "<<fq<<" "<<qtotal<<" "<<qavg<<" "<<qclus<<" "<<qscale<<" "
3137  // <<fbin_[0]<<" "<<fbin_[1]<<" "<<fbin_[2]<<std::endl;
3138  } else {
3139  if (fq > fbin_[1]) {
3140  binq = 1;
3141  } else {
3142  if (fq > fbin_[2]) {
3143  binq = 2;
3144  } else {
3145  binq = 3;
3146  }
3147  }
3148  }
3149 
3150  auto yavggen = (1.f - yratio_) * enty0_->yavggen[binq] + yratio_ * enty1_->yavggen[binq];
3151  if (flip_y_) {
3152  yavggen = -yavggen;
3153  }
3154  auto yrmsgen = (1.f - yratio_) * enty0_->yrmsgen[binq] + yratio_ * enty1_->yrmsgen[binq];
3155 
3156  // next, loop over all x-angle entries, first, find relevant y-slices
3157 
3158  auto iylow = 0;
3159  auto iyhigh = 0;
3160  auto yxratio = 0.f;
3161 
3162  {
3163  auto j = std::lower_bound(templ.cotbetaX.begin(), templ.cotbetaX.begin() + Nyx, acotb);
3164  if (j == templ.cotbetaX.begin() + Nyx) {
3165  --j;
3166  yxratio = 1.f;
3167  } else if (j == templ.cotbetaX.begin()) {
3168  ++j;
3169  yxratio = 0.f;
3170  } else {
3171  yxratio = (acotb - (*(j - 1))) / ((*j) - (*(j - 1)));
3172  }
3173 
3174  iyhigh = j - templ.cotbetaX.begin();
3175  iylow = iyhigh - 1;
3176  }
3177 
3178  auto ilow = 0;
3179  auto ihigh = 0;
3180  auto xxratio = 0.f;
3181 
3182  {
3183  auto j = std::lower_bound(templ.cotalphaX.begin(), templ.cotalphaX.begin() + Nxx, cota);
3184  if (j == templ.cotalphaX.begin() + Nxx) {
3185  --j;
3186  xxratio = 1.f;
3187  } else if (j == templ.cotalphaX.begin()) {
3188  ++j;
3189  xxratio = 0.f;
3190  } else {
3191  xxratio = (cota - (*(j - 1))) / ((*j) - (*(j - 1)));
3192  }
3193 
3194  ihigh = j - templ.cotalphaX.begin();
3195  ilow = ihigh - 1;
3196  }
3197 
3198  dx1 =
3199  (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].dxone + xxratio * thePixelTemp_[index].entx[0][ihigh].dxone;
3200  if (flip_x) {
3201  dx1 = -dx1;
3202  }
3203  sx1 =
3204  (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].sxone + xxratio * thePixelTemp_[index].entx[0][ihigh].sxone;
3205  dx2 =
3206  (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].dxtwo + xxratio * thePixelTemp_[index].entx[0][ihigh].dxtwo;
3207  if (flip_x) {
3208  dx2 = -dx2;
3209  }
3210  sx2 =
3211  (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].sxtwo + xxratio * thePixelTemp_[index].entx[0][ihigh].sxtwo;
3212 
3213  // pixmax is the maximum allowed pixel charge (used for truncation)
3214 
3215  pixmx = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].pixmax +
3216  xxratio * thePixelTemp_[index].entx[iylow][ihigh].pixmax) +
3217  yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].pixmax +
3218  xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].pixmax);
3219 
3220  auto xavggen = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].xavggen[binq] +
3221  xxratio * thePixelTemp_[index].entx[iylow][ihigh].xavggen[binq]) +
3222  yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].xavggen[binq] +
3223  xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].xavggen[binq]);
3224  if (flip_x) {
3225  xavggen = -xavggen;
3226  }
3227 
3228  auto xrmsgen = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].xrmsgen[binq] +
3229  xxratio * thePixelTemp_[index].entx[iylow][ihigh].xrmsgen[binq]) +
3230  yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].xrmsgen[binq] +
3231  xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].xrmsgen[binq]);
3232 
3233  // Take the errors and bias from the correct charge bin
3234 
3235  sigmay = yrmsgen;
3236  deltay = yavggen;
3237 
3238  sigmax = xrmsgen;
3239  deltax = xavggen;
3240 
3241  // If the charge is too small (then flag it)
3242 
3243  if (qtotal < 0.95f * qmin) {
3244  binq = 5;
3245  } else {
3246  if (qtotal < 0.95f * qmin2) {
3247  binq = 4;
3248  }
3249  }
3250 
3251  return binq;
3252 
3253 } // qbin
3254 
3255 // ************************************************************************************************************
3283 // ************************************************************************************************************
3284 int SiPixelTemplate::qbin(int id,
3285  float cotalpha,
3286  float cotbeta,
3287  float locBz,
3288  float qclus,
3289  float& pixmx,
3290  float& sigmay,
3291  float& deltay,
3292  float& sigmax,
3293  float& deltax,
3294  float& sy1,
3295  float& dy1,
3296  float& sy2,
3297  float& dy2,
3298  float& sx1,
3299  float& dx1,
3300  float& sx2,
3301  float& dx2) // float& lorywidth, float& lorxwidth)
3302 {
3303  // Interpolate for a new set of track angles
3304 
3305  // Local variables
3306  float locBx = 1.f; // lorywidth, lorxwidth;
3307 
3308  return SiPixelTemplate::qbin(id,
3309  cotalpha,
3310  cotbeta,
3311  locBz,
3312  locBx,
3313  qclus,
3314  pixmx,
3315  sigmay,
3316  deltay,
3317  sigmax,
3318  deltax,
3319  sy1,
3320  dy1,
3321  sy2,
3322  dy2,
3323  sx1,
3324  dx1,
3325  sx2,
3326  dx2); // , lorywidth, lorxwidth);
3327 
3328 } // qbin
3329 
3330 // ************************************************************************************************************
3352 // ************************************************************************************************************
3353 /*
3354  int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float& pixmx, float& sigmay, float& deltay, float& sigmax, float& deltax,
3355  float& sy1, float& dy1, float& sy2, float& dy2, float& sx1, float& dx1, float& sx2, float& dx2)
3356 
3357  {
3358  float lorywidth, lorxwidth;
3359  return SiPixelTemplate::qbin(id, cotalpha, cotbeta, locBz, qclus, pixmx, sigmay, deltay, sigmax, deltax,
3360  sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, lorywidth, lorxwidth);
3361 
3362  } // qbin
3363  */
3364 
3365 // ************************************************************************************************************
3372 // ************************************************************************************************************
3373 int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float qclus) {
3374  // Interpolate for a new set of track angles
3375 
3376  // Local variables
3377  float pixmx, sigmay, deltay, sigmax, deltax, sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, locBz; // lorywidth, lorxwidth;
3378  // Local variables
3379  float locBx = 1.f;
3380  if (cotbeta < 0.f) {
3381  locBx = -1.f;
3382  }
3383  locBz = locBx;
3384  if (cotalpha < 0.f) {
3385  locBz = -locBx;
3386  }
3387 
3388  return SiPixelTemplate::qbin(id,
3389  cotalpha,
3390  cotbeta,
3391  locBz,
3392  locBx,
3393  qclus,
3394  pixmx,
3395  sigmay,
3396  deltay,
3397  sigmax,
3398  deltax,
3399  sy1,
3400  dy1,
3401  sy2,
3402  dy2,
3403  sx1,
3404  dx1,
3405  sx2,
3406  dx2); // , lorywidth, lorxwidth);
3407 
3408 } // qbin
3409 
3410 // ************************************************************************************************************
3416 // ************************************************************************************************************
3417 int SiPixelTemplate::qbin(int id, float cotbeta, float qclus) {
3418  // Interpolate for a new set of track angles
3419 
3420  // Local variables
3421  float pixmx, sigmay, deltay, sigmax, deltax, sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, locBz,
3422  locBx; //, lorywidth, lorxwidth;
3423  const float cotalpha = 0.f;
3424  locBx = 1.f;
3425  if (cotbeta < 0.f) {
3426  locBx = -1.f;
3427  }
3428  locBz = locBx;
3429  if (cotalpha < 0.f) {
3430  locBz = -locBx;
3431  }
3432  return SiPixelTemplate::qbin(id,
3433  cotalpha,
3434  cotbeta,
3435  locBz,
3436  locBx,
3437  qclus,
3438  pixmx,
3439  sigmay,
3440  deltay,
3441  sigmax,
3442  deltax,
3443  sy1,
3444  dy1,
3445  sy2,
3446  dy2,
3447  sx1,
3448  dx1,
3449  sx2,
3450  dx2); // , lorywidth, lorxwidth);
3451 
3452 } // qbin
3453 
3454 // ************************************************************************************************************
3466 // ************************************************************************************************************
3467 void SiPixelTemplate::temperrors(int id,
3468  float cotalpha,
3469  float cotbeta,
3470  int qBin,
3471  float& sigmay,
3472  float& sigmax,
3473  float& sy1,
3474  float& sy2,
3475  float& sx1,
3476  float& sx2)
3477 
3478 {
3479  // Interpolate for a new set of track angles
3480 
3481  // Local variables
3482  int i;
3483  int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, index;
3484  float yxratio, xxratio;
3485  float acotb;
3486  float yrms, xrms;
3487  //bool flip_y;
3488 
3489  // Find the index corresponding to id
3490 
3491  index = -1;
3492  for (i = 0; i < (int)thePixelTemp_.size(); ++i) {
3493  if (id == thePixelTemp_[i].head.ID) {
3494  index = i;
3495  break;
3496  }
3497  }
3498 
3499 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3500  if (index < 0 || index >= (int)thePixelTemp_.size()) {
3501  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors can't find needed template ID = " << id
3502  << std::endl;
3503  }
3504 #else
3505  assert(index >= 0 && index < (int)thePixelTemp_.size());
3506 #endif
3507 
3508 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3509  if (qBin < 0 || qBin > 5) {
3510  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors called with illegal qBin = " << qBin
3511  << std::endl;
3512  }
3513 #else
3514  assert(qBin >= 0 && qBin < 6);
3515 #endif
3516 
3517  // The error information for qBin > 3 is taken to be the same as qBin=3
3518 
3519  if (qBin > 3) {
3520  qBin = 3;
3521  }
3522  //
3523 
3524  // Interpolate the absolute value of cot(beta)
3525 
3526  acotb = std::abs(cotbeta);
3527 
3528  // Copy the charge scaling factor to the private variable
3529 
3530  Ny = thePixelTemp_[index].head.NTy;
3531  Nyx = thePixelTemp_[index].head.NTyx;
3532  Nxx = thePixelTemp_[index].head.NTxx;
3533 
3534 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3535  if (Ny < 2 || Nyx < 1 || Nxx < 2) {
3536  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny
3537  << "/" << Nyx << "/" << Nxx << std::endl;
3538  }
3539 #else
3540  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
3541 #endif
3542 
3543  // next, loop over all y-angle entries
3544 
3545  // Interpolate/store all y-related quantities (flip displacements when flip_y)
3546 
3547  sy1 = (1.f - yratio_) * enty0_->syone + yratio_ * enty1_->syone;
3548  sy2 = (1.f - yratio_) * enty0_->sytwo + yratio_ * enty1_->sytwo;
3549  yrms = (1.f - yratio_) * enty0_->yrms[qBin] + yratio_ * enty1_->yrms[qBin];
3550 
3551  // next, loop over all x-angle entries, first, find relevant y-slices
3552 
3553  iylow = 0;
3554  yxratio = 0.f;
3555 
3556  if (acotb >= thePixelTemp_[index].entx[Nyx - 1][0].cotbeta) {
3557  iylow = Nyx - 2;
3558  yxratio = 1.f;
3559 
3560  } else if (acotb >= thePixelTemp_[index].entx[0][0].cotbeta) {
3561  for (i = 0; i < Nyx - 1; ++i) {
3562  if (thePixelTemp_[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp_[index].entx[i + 1][0].cotbeta) {
3563  iylow = i;
3564  yxratio = (acotb - thePixelTemp_[index].entx[i][0].cotbeta) /
3565  (thePixelTemp_[index].entx[i + 1][0].cotbeta - thePixelTemp_[index].entx[i][0].cotbeta);
3566  break;
3567  }
3568  }
3569  }
3570 
3571  iyhigh = iylow + 1;
3572 
3573  ilow = 0;
3574  xxratio = 0.f;
3575 
3576  if (cotalpha >= thePixelTemp_[index].entx[0][Nxx - 1].cotalpha) {
3577  ilow = Nxx - 2;
3578  xxratio = 1.f;
3579 
3580  } else {
3581  if (cotalpha >= thePixelTemp_[index].entx[0][0].cotalpha) {
3582  for (i = 0; i < Nxx - 1; ++i) {
3583  if (thePixelTemp_[index].entx[0][i].cotalpha <= cotalpha &&
3584  cotalpha < thePixelTemp_[index].entx[0][i + 1].cotalpha) {
3585  ilow = i;
3586  xxratio = (cotalpha - thePixelTemp_[index].entx[0][i].cotalpha) /
3587  (thePixelTemp_[index].entx[0][i + 1].cotalpha - thePixelTemp_[index].entx[0][i].cotalpha);
3588  break;
3589  }
3590  }
3591  }
3592  }
3593 
3594  ihigh = ilow + 1;
3595 
3596  sx1 =
3597  (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].sxone + xxratio * thePixelTemp_[index].entx[0][ihigh].sxone;
3598  sx2 =
3599  (1.f - xxratio) * thePixelTemp_[index].entx[0][ilow].sxtwo + xxratio * thePixelTemp_[index].entx[0][ihigh].sxtwo;
3600 
3601  xrms = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].xrms[qBin] +
3602  xxratio * thePixelTemp_[index].entx[iylow][ihigh].xrms[qBin]) +
3603  yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].xrms[qBin] +
3604  xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].xrms[qBin]);
3605 
3606  // Take the errors and bias from the correct charge bin
3607 
3608  sigmay = yrms;
3609 
3610  sigmax = xrms;
3611 
3612  return;
3613 
3614 } // temperrors
3615 
3616 // ************************************************************************************************************
3626 // ************************************************************************************************************
3627 void SiPixelTemplate::qbin_dist(int id,
3628  float cotalpha,
3629  float cotbeta,
3630  float qbin_frac[4],
3631  float& ny1_frac,
3632  float& ny2_frac,
3633  float& nx1_frac,
3634  float& nx2_frac)
3635 
3636 {
3637  // Interpolate for a new set of track angles
3638 
3639  // Local variables
3640  int i;
3641  int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, index;
3642  float yxratio, xxratio;
3643  float acotb;
3644  float qfrac[4];
3645  //bool flip_y;
3646 
3647  // Find the index corresponding to id
3648 
3649  index = -1;
3650  for (i = 0; i < (int)thePixelTemp_.size(); ++i) {
3651  if (id == thePixelTemp_[i].head.ID) {
3652  index = i;
3653  // id_current_ = id;
3654  break;
3655  }
3656  }
3657 
3658 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3659  if (index < 0 || index >= (int)thePixelTemp_.size()) {
3660  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors can't find needed template ID = " << id
3661  << std::endl;
3662  }
3663 #else
3664  assert(index >= 0 && index < (int)thePixelTemp_.size());
3665 #endif
3666 
3667  //
3668 
3669  // Interpolate the absolute value of cot(beta)
3670 
3671  acotb = fabs((double)cotbeta);
3672 
3673  // Copy the charge scaling factor to the private variable
3674 
3675  Ny = thePixelTemp_[index].head.NTy;
3676  Nyx = thePixelTemp_[index].head.NTyx;
3677  Nxx = thePixelTemp_[index].head.NTxx;
3678 
3679 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3680  if (Ny < 2 || Nyx < 1 || Nxx < 2) {
3681  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny
3682  << "/" << Nyx << "/" << Nxx << std::endl;
3683  }
3684 #else
3685  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
3686 #endif
3687 
3688  // Interpolate/store all y-related quantities (flip displacements when flip_y)
3689  ny1_frac = (1.f - yratio_) * enty0_->fracyone + yratio_ * enty1_->fracyone;
3690  ny2_frac = (1.f - yratio_) * enty0_->fracytwo + yratio_ * enty1_->fracytwo;
3691 
3692  // next, loop over all x-angle entries, first, find relevant y-slices
3693 
3694  iylow = 0;
3695  yxratio = 0.f;
3696 
3697  if (acotb >= thePixelTemp_[index].entx[Nyx - 1][0].cotbeta) {
3698  iylow = Nyx - 2;
3699  yxratio = 1.f;
3700 
3701  } else if (acotb >= thePixelTemp_[index].entx[0][0].cotbeta) {
3702  for (i = 0; i < Nyx - 1; ++i) {
3703  if (thePixelTemp_[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp_[index].entx[i + 1][0].cotbeta) {
3704  iylow = i;
3705  yxratio = (acotb - thePixelTemp_[index].entx[i][0].cotbeta) /
3706  (thePixelTemp_[index].entx[i + 1][0].cotbeta - thePixelTemp_[index].entx[i][0].cotbeta);
3707  break;
3708  }
3709  }
3710  }
3711 
3712  iyhigh = iylow + 1;
3713 
3714  ilow = 0;
3715  xxratio = 0.f;
3716 
3717  if (cotalpha >= thePixelTemp_[index].entx[0][Nxx - 1].cotalpha) {
3718  ilow = Nxx - 2;
3719  xxratio = 1.f;
3720 
3721  } else {
3722  if (cotalpha >= thePixelTemp_[index].entx[0][0].cotalpha) {
3723  for (i = 0; i < Nxx - 1; ++i) {
3724  if (thePixelTemp_[index].entx[0][i].cotalpha <= cotalpha &&
3725  cotalpha < thePixelTemp_[index].entx[0][i + 1].cotalpha) {
3726  ilow = i;
3727  xxratio = (cotalpha - thePixelTemp_[index].entx[0][i].cotalpha) /
3728  (thePixelTemp_[index].entx[0][i + 1].cotalpha - thePixelTemp_[index].entx[0][i].cotalpha);
3729  break;
3730  }
3731  }
3732  }
3733  }
3734 
3735  ihigh = ilow + 1;
3736 
3737  for (i = 0; i < 3; ++i) {
3738  qfrac[i] = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].qbfrac[i] +
3739  xxratio * thePixelTemp_[index].entx[iylow][ihigh].qbfrac[i]) +
3740  yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].qbfrac[i] +
3741  xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].qbfrac[i]);
3742  }
3743  nx1_frac = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].fracxone +
3744  xxratio * thePixelTemp_[index].entx[iylow][ihigh].fracxone) +
3745  yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].fracxone +
3746  xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].fracxone);
3747  nx2_frac = (1.f - yxratio) * ((1.f - xxratio) * thePixelTemp_[index].entx[iylow][ilow].fracxtwo +
3748  xxratio * thePixelTemp_[index].entx[iylow][ihigh].fracxtwo) +
3749  yxratio * ((1.f - xxratio) * thePixelTemp_[index].entx[iyhigh][ilow].fracxtwo +
3750  xxratio * thePixelTemp_[index].entx[iyhigh][ihigh].fracxtwo);
3751 
3752  qbin_frac[0] = qfrac[0];
3753  qbin_frac[1] = qbin_frac[0] + qfrac[1];
3754  qbin_frac[2] = qbin_frac[1] + qfrac[2];
3755  qbin_frac[3] = 1.f;
3756  return;
3757 
3758 } // qbin_dist
3759 
3760 // *************************************************************************************************************************************
3762 
3768 // *************************************************************************************************************************************
3769 
3770 bool SiPixelTemplate::simpletemplate2D(
3771  float xhit, float yhit, std::vector<bool>& ydouble, std::vector<bool>& xdouble, float template2d[BXM2][BYM2]) {
3772  // Local variables
3773 
3774  float x0, y0, xf, yf, xi, yi, sf, si, s0, qpix, slopey, slopex, ds;
3775  int i, j, jpix0, ipix0, jpixf, ipixf, jpix, ipix, nx, ny, anx, any, jmax, imax;
3776  float qtotal;
3777  // double path;
3778  std::list<SimplePixel> list;
3779  std::list<SimplePixel>::iterator listIter, listEnd;
3780 
3781  // Calculate the entry and exit points for the line charge from the track
3782 
3783  x0 = xhit - 0.5 * zsize_ * cota_current_;
3784  y0 = yhit - 0.5 * zsize_ * cotb_current_;
3785 
3786  jpix0 = floor(x0 / xsize_) + 1;
3787  ipix0 = floor(y0 / ysize_) + 1;
3788 
3789  if (jpix0 < 0 || jpix0 > BXM3) {
3790  return false;
3791  }
3792  if (ipix0 < 0 || ipix0 > BYM3) {
3793  return false;
3794  }
3795 
3796  xf = xhit + 0.5 * zsize_ * cota_current_ + lorxwidth_;
3797  yf = yhit + 0.5 * zsize_ * cotb_current_ + lorywidth_;
3798 
3799  jpixf = floor(xf / xsize_) + 1;
3800  ipixf = floor(yf / ysize_) + 1;
3801 
3802  if (jpixf < 0 || jpixf > BXM3) {
3803  return false;
3804  }
3805  if (ipixf < 0 || ipixf > BYM3) {
3806  return false;
3807  }
3808 
3809  // total charge length
3810 
3811  sf = std::sqrt((xf - x0) * (xf - x0) + (yf - y0) * (yf - y0));
3812  if ((xf - x0) != 0.f) {
3813  slopey = (yf - y0) / (xf - x0);
3814  } else {
3815  slopey = 1.e10;
3816  }
3817  if ((yf - y0) != 0.f) {
3818  slopex = (xf - x0) / (yf - y0);
3819  } else {
3820  slopex = 1.e10;
3821  }
3822 
3823  // use average charge in this direction
3824 
3825  qtotal = qavg_avg_;
3826 
3827  SimplePixel element;
3828  element.s = sf;
3829  element.x = xf;
3830  element.y = yf;
3831  element.i = ipixf;
3832  element.j = jpixf;
3833  element.btype = 0;
3834  list.push_back(element);
3835 
3836  // nx is the number of x interfaces crossed by the line charge
3837 
3838  nx = jpixf - jpix0;
3839  anx = abs(nx);
3840  if (anx > 0) {
3841  if (nx > 0) {
3842  for (j = jpix0; j < jpixf; ++j) {
3843  xi = xsize_ * j;
3844  yi = slopey * (xi - x0) + y0;
3845  ipix = (int)(yi / ysize_) + 1;
3846  si = std::sqrt((xi - x0) * (xi - x0) + (yi - y0) * (yi - y0));
3847  element.s = si;
3848  element.x = xi;
3849  element.y = yi;
3850  element.i = ipix;
3851  element.j = j;
3852  element.btype = 1;
3853  list.push_back(element);
3854  }
3855  } else {
3856  for (j = jpix0; j > jpixf; --j) {
3857  xi = xsize_ * (j - 1);
3858  yi = slopey * (xi - x0) + y0;
3859  ipix = (int)(yi / ysize_) + 1;
3860  si = std::sqrt((xi - x0) * (xi - x0) + (yi - y0) * (yi - y0));
3861  element.s = si;
3862  element.x = xi;
3863  element.y = yi;
3864  element.i = ipix;
3865  element.j = j;
3866  element.btype = 1;
3867  list.push_back(element);
3868  }
3869  }
3870  }
3871 
3872  ny = ipixf - ipix0;
3873  any = abs(ny);
3874  if (any > 0) {
3875  if (ny > 0) {
3876  for (i = ipix0; i < ipixf; ++i) {
3877  yi = ysize_ * i;
3878  xi = slopex * (yi - y0) + x0;
3879  jpix = (int)(xi / xsize_) + 1;
3880  si = std::sqrt((xi - x0) * (xi - x0) + (yi - y0) * (yi - y0));
3881  element.s = si;
3882  element.x = xi;
3883  element.y = yi;
3884  element.i = i;
3885  element.j = jpix;
3886  element.btype = 2;
3887  list.push_back(element);
3888  }
3889  } else {
3890  for (i = ipix0; i > ipixf; --i) {
3891  yi = ysize_ * (i - 1);
3892  xi = slopex * (yi - y0) + x0;
3893  jpix = (int)(xi / xsize_) + 1;
3894  si = std::sqrt((xi - x0) * (xi - x0) + (yi - y0) * (yi - y0));
3895  element.s = si;
3896  element.x = xi;
3897  element.y = yi;
3898  element.i = i;
3899  element.j = jpix;
3900  element.btype = 2;
3901  list.push_back(element);
3902  }
3903  }
3904  }
3905 
3906  imax = std::max(ipix0, ipixf);
3907  jmax = std::max(jpix0, jpixf);
3908 
3909  // Sort the list according to the distance from the initial point
3910 
3911  list.sort();
3912 
3913  // Look for double pixels and adjust the list appropriately
3914 
3915  for (i = 1; i < imax; ++i) {
3916  if (ydouble[i - 1]) {
3917  listIter = list.begin();
3918  if (ny > 0) {
3919  while (listIter != list.end()) {
3920  if (listIter->i == i && listIter->btype == 2) {
3921  listIter = list.erase(listIter);
3922  continue;
3923  }
3924  if (listIter->i > i) {
3925  --(listIter->i);
3926  }
3927  ++listIter;
3928  }
3929  } else {
3930  while (listIter != list.end()) {
3931  if (listIter->i == i + 1 && listIter->btype == 2) {
3932  listIter = list.erase(listIter);
3933  continue;
3934  }
3935  if (listIter->i > i + 1) {
3936  --(listIter->i);
3937  }
3938  ++listIter;
3939  }
3940  }
3941  }
3942  }
3943 
3944  for (j = 1; j < jmax; ++j) {
3945  if (xdouble[j - 1]) {
3946  listIter = list.begin();
3947  if (nx > 0) {
3948  while (listIter != list.end()) {
3949  if (listIter->j == j && listIter->btype == 1) {
3950  listIter = list.erase(listIter);
3951  continue;
3952  }
3953  if (listIter->j > j) {
3954  --(listIter->j);
3955  }
3956  ++listIter;
3957  }
3958  } else {
3959  while (listIter != list.end()) {
3960  if (listIter->j == j + 1 && listIter->btype == 1) {
3961  listIter = list.erase(listIter);
3962  continue;
3963  }
3964  if (listIter->j > j + 1) {
3965  --(listIter->j);
3966  }
3967  ++listIter;
3968  }
3969  }
3970  }
3971  }
3972 
3973  // The list now contains the path lengths of the line charge in each pixel from (x0,y0). Cacluate the lengths of the segments and the charge.
3974 
3975  s0 = 0.f;
3976  listIter = list.begin();
3977  listEnd = list.end();
3978  for (; listIter != listEnd; ++listIter) {
3979  si = listIter->s;
3980  ds = si - s0;
3981  s0 = si;
3982  j = listIter->j;
3983  i = listIter->i;
3984  if (sf > 0.f) {
3985  qpix = qtotal * ds / sf;
3986  } else {
3987  qpix = qtotal;
3988  }
3989  template2d[j][i] += qpix;
3990  }
3991 
3992  return true;
3993 
3994 } // simpletemplate2D
3995 
3996 // ************************************************************************************************************
4001 // ************************************************************************************************************
4002 void SiPixelTemplate::vavilov_pars(double& mpv, double& sigma, double& kappa)
4003 
4004 {
4005  // Local variables
4006  int i;
4007  int ilow, ihigh, Ny;
4008  float yratio, cotb, cotalpha0, arg;
4009 
4010  // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta)
4011 
4012  cotalpha0 = thePixelTemp_[index_id_].enty[0].cotalpha;
4013  arg = cotb_current_ * cotb_current_ + cota_current_ * cota_current_ - cotalpha0 * cotalpha0;
4014  if (arg < 0.f)
4015  arg = 0.f;
4016  cotb = std::sqrt(arg);
4017 
4018  // Copy the charge scaling factor to the private variable
4019 
4020  Ny = thePixelTemp_[index_id_].head.NTy;
4021 
4022 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
4023  if (Ny < 2) {
4024  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny
4025  << std::endl;
4026  }
4027 #else
4028  assert(Ny > 1);
4029 #endif
4030 
4031  // next, loop over all y-angle entries
4032 
4033  ilow = 0;
4034  yratio = 0.f;
4035 
4036  if (cotb >= thePixelTemp_[index_id_].enty[Ny - 1].cotbeta) {
4037  ilow = Ny - 2;
4038  yratio = 1.f;
4039 
4040  } else {
4041  if (cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
4042  for (i = 0; i < Ny - 1; ++i) {
4043  if (thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i + 1].cotbeta) {
4044  ilow = i;
4045  yratio = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta) /
4046  (thePixelTemp_[index_id_].enty[i + 1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
4047  break;
4048  }
4049  }
4050  }
4051  }
4052 
4053  ihigh = ilow + 1;
4054 
4055  // Interpolate Vavilov parameters
4056 
4057  mpvvav_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].mpvvav +
4058  yratio * thePixelTemp_[index_id_].enty[ihigh].mpvvav;
4059  sigmavav_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].sigmavav +
4060  yratio * thePixelTemp_[index_id_].enty[ihigh].sigmavav;
4061  kappavav_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].kappavav +
4062  yratio * thePixelTemp_[index_id_].enty[ihigh].kappavav;
4063 
4064  // Copy to parameter list
4065 
4066  //Avoid rounding difference between floats and doubles causing issues later
4067  if (kappavav_ <= 0.01f) {
4068  LOGERROR("SiPixelTemplate") << "Vavilov kappa value is " << kappavav_ << " changing it to be above 0.01" << ENDL;
4069  kappavav_ = 0.01f + 0.0000001f;
4070  }
4071 
4072  mpv = (double)mpvvav_;
4073  sigma = (double)sigmavav_;
4074  kappa = (double)kappavav_;
4075 
4076  return;
4077 
4078 } // vavilov_pars
4079 
4080 // ************************************************************************************************************
4085 // ************************************************************************************************************
4086 void SiPixelTemplate::vavilov2_pars(double& mpv, double& sigma, double& kappa)
4087 
4088 {
4089  // Local variables
4090  int i;
4091  int ilow, ihigh, Ny;
4092  float yratio, cotb, cotalpha0, arg;
4093 
4094  // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta)
4095 
4096  cotalpha0 = thePixelTemp_[index_id_].enty[0].cotalpha;
4097  arg = cotb_current_ * cotb_current_ + cota_current_ * cota_current_ - cotalpha0 * cotalpha0;
4098  if (arg < 0.f)
4099  arg = 0.f;
4100  cotb = std::sqrt(arg);
4101 
4102  // Copy the charge scaling factor to the private variable
4103 
4104  Ny = thePixelTemp_[index_id_].head.NTy;
4105 
4106 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
4107  if (Ny < 2) {
4108  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny
4109  << std::endl;
4110  }
4111 #else
4112  assert(Ny > 1);
4113 #endif
4114 
4115  // next, loop over all y-angle entries
4116 
4117  ilow = 0;
4118  yratio = 0.f;
4119 
4120  if (cotb >= thePixelTemp_[index_id_].enty[Ny - 1].cotbeta) {
4121  ilow = Ny - 2;
4122  yratio = 1.f;
4123 
4124  } else {
4125  if (cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
4126  for (i = 0; i < Ny - 1; ++i) {
4127  if (thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i + 1].cotbeta) {
4128  ilow = i;
4129  yratio = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta) /
4130  (thePixelTemp_[index_id_].enty[i + 1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
4131  break;
4132  }
4133  }
4134  }
4135  }
4136 
4137  ihigh = ilow + 1;
4138 
4139  // Interpolate Vavilov parameters
4140 
4141  mpvvav2_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].mpvvav2 +
4142  yratio * thePixelTemp_[index_id_].enty[ihigh].mpvvav2;
4143  sigmavav2_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].sigmavav2 +
4144  yratio * thePixelTemp_[index_id_].enty[ihigh].sigmavav2;
4145  kappavav2_ = (1.f - yratio) * thePixelTemp_[index_id_].enty[ilow].kappavav2 +
4146  yratio * thePixelTemp_[index_id_].enty[ihigh].kappavav2;
4147 
4148  // Copy to parameter list
4149 
4150  mpv = (double)mpvvav2_;
4151  sigma = (double)sigmavav2_;
4152  kappa = (double)kappavav2_;
4153 
4154  return;
4155 
4156 } // vavilov2_pars
SiPixelTemplate::cota_current_
float cota_current_
current cot alpha
Definition: SiPixelTemplate.h:691
SiPixelTemplateEntry::qavg_avg
float qavg_avg
average cluster charge of clusters that are less than qavg (normalize 2-D simple templates) ...
Definition: SiPixelTemplate.h:177
SiPixelTemplate::mpvvav_
float mpvvav_
most probable charge in Vavilov distribution (not actually for larger kappa)
Definition: SiPixelTemplate.h:773
SiPixelTemplateEntry::runnum
int runnum
&lt; Basic template entry corresponding to a single set of track angles
Definition: SiPixelTemplate.h:105
SiPixelTemplateEntry::xflpar
float xflpar[4][6]
Aqfl-parameterized x-correction in 4 charge bins.
Definition: SiPixelTemplate.h:151
SiPixelTemplateEntry::xrms
float xrms[4]
average x-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:148
SiPixelTemplate::vavilov2_pars
void vavilov2_pars(double &mpv, double &sigma, double &kappa)
Definition: SiPixelTemplate.cc:4086
SiPixelTemplateStore::cotbetaY
std::array< float, 100 > cotbetaY
Definition: SiPixelTemplate.h:210
BXSIZE
#define BXSIZE
Definition: SiPixelTemplateDefs.h:33
LOGERROR
#define LOGERROR(x)
Definition: SiPixelTemplate.cc:110
SiPixelTemplateEntry::clslenx
float clslenx
cluster x-length in pixels at signal height sxmax/2
Definition: SiPixelTemplate.h:131
c
const edm::EventSetup & c
Definition: SiStripLAProfileBooker.cc:66
SiPixelTemplateEntry::qavg
float qavg
average cluster charge for this set of track angles (now includes threshold effects) ...
Definition: SiPixelTemplate.h:111
mps_fire.i
i
Definition: mps_fire.py:428
SiPixelTemplateEntry::syone
float syone
rms for one pixel y-clusters
Definition: SiPixelTemplate.h:115
SiPixelTemplateEntry::chi2yavgone
float chi2yavgone
average y chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:156
SiPixelTemplate::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
Definition: SiPixelTemplate.cc:1328
SiPixelTemplate::kappavav2
float kappavav2()
kappa parameter for 2-cluster Vavilov distribution
Definition: SiPixelTemplate.h:674
SiPixelTemplate::cytemp
int cytemp()
Return central pixel of y template pixels above readout threshold.
Definition: SiPixelTemplate.cc:2656
gpuClustering::id
uint16_t *__restrict__ id
Definition: gpuClusterChargeCut.h:21
SiPixelTemplateEntry::dyone
float dyone
mean offset/correction for one pixel y-clusters
Definition: SiPixelTemplate.h:114
BYSIZE
#define BYSIZE
Definition: SiPixelTemplateDefs.h:25
SiPixelTemplateEntry::sigmavav
float sigmavav
&quot;sigma&quot; scale fctor for Vavilov distribution
Definition: SiPixelTemplate.h:133
SiPixelTemplateEntry::fracxtwo
float fracxtwo
fraction of double pixel sample with xsize = 1
Definition: SiPixelTemplate.h:176
SiPixelTemplateEntry::yavggen
float yavggen[4]
generic algorithm: average y-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:125
SiPixelTemplateEntry::sxmax
float sxmax
average pixel signal for x-projection of cluster
Definition: SiPixelTemplate.h:116
SimplePixel::btype
int btype
type of boundary (0=end, 1 = x-boundary, 2 = y-boundary)
Definition: SimplePixel.h:29
SiPixelTemplate::xsigma2
void xsigma2(int fxpix, int lxpix, float sxthr, float xsum[13+4], float xsig2[13+4])
Definition: SiPixelTemplate.cc:2274
SiPixelTemplateEntry::chi2yminc2m
float chi2yminc2m[4]
1st pass chi2 min search: minimum of y chi^2 in 4 charge bins (merged clusters)
Definition: SiPixelTemplate.h:163
SiPixelTemplateEntry::xavgc2m
float xavgc2m[4]
1st pass chi2 min search: average x-bias of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:164
TXSIZE
#define TXSIZE
Definition: SiPixelTemplateDefs.h:30
any
bool any(const std::vector< T > &v, const T &what)
Definition: ECalSD.cc:37
SiPixelTemplateEntry::chi2xavgone
float chi2xavgone
average x chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:158
SiPixelTemplateEntry::fracytwo
float fracytwo
fraction of double pixel sample with ysize = 1
Definition: SiPixelTemplate.h:175
SiPixelTemplateEntry::yrms
float yrms[4]
average y-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:143
SiPixelTemplateStore::cotalphaX
std::array< float, 80 > cotalphaX
60 y templates spanning cluster lengths from 0px to +18px
Definition: SiPixelTemplate.h:212
SiPixelTemplateHeader::NTxx
int NTxx
number of Template x-entries in each slice
Definition: SiPixelTemplate.h:186
BXM3
#define BXM3
Definition: SiPixelTemplateDefs.h:37
BXM1
#define BXM1
Definition: SiPixelTemplateDefs.h:35
SiPixelTemplate::mpvvav
float mpvvav()
most probable charge in Vavilov distribution (not actually for larger kappa)
Definition: SiPixelTemplate.h:667
SiPixelTemplateEntry::mpvvav2
float mpvvav2
most probable charge in Vavilov distribution for 2 merged clusters (not actually for larger kappa) ...
Definition: SiPixelTemplate.h:135
SiPixelTemplateEntry::ygx0gen
float ygx0gen[4]
generic algorithm: average y0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:168
SiPixelTemplateEntry::xavg
float xavg[4]
average x-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:147
EcalCondDB.db
tuple db
Definition: EcalCondDB.py:153
SiPixelTemplate::simpletemplate2D
bool simpletemplate2D(float xhitp, float yhitp, std::vector< bool > &ydouble, std::vector< bool > &xdouble, float template2d[13+2][21+2])
Make simple 2-D templates from track angles set in interpolate and hit position.
Definition: SiPixelTemplate.cc:3770
SiPixelTemplateHeader::lorybias
float lorybias
estimate of y-lorentz bias
Definition: SiPixelTemplate.h:191
SimplePixel::y
float y
y coordinate of boundary intersection
Definition: SimplePixel.h:26
SiPixelTemplate::yratio
float yratio()
fractional distance in y between cotbeta templates
Definition: SiPixelTemplate.h:439
SiPixelTemplateHeader::fluence
float fluence
radiation fluence in n_eq/cm^2
Definition: SiPixelTemplate.h:195
SiPixelTemplateHeader::ss50
float ss50
1/2 of the single hit dcol threshold in electrons
Definition: SiPixelTemplate.h:197
SiPixelTemplateEntry::dytwo
float dytwo
mean offset/correction for one double-pixel y-clusters
Definition: SiPixelTemplate.h:119
cond::impl::to_string
std::string to_string(const V &value)
Definition: OMSAccess.h:71
SiPixelTemplate::xrms
float xrms(int i)
average x-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:492
SiPixelTemplateHeader::ysize
float ysize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate.h:203
SiPixelTemplateEntry::dxone
float dxone
mean offset/correction for one pixel x-clusters
Definition: SiPixelTemplate.h:117
SiPixelTemplateEntry::pixmax
float pixmax
maximum charge for individual pixels in cluster
Definition: SiPixelTemplate.h:112
SiPixelTemplateHeader::qscale
float qscale
Charge scaling to match cmssw and pixelav.
Definition: SiPixelTemplate.h:188
SiPixelTemplateEntry::chi2yavgc2m
float chi2yavgc2m[4]
1st pass chi2 min search: average y chi^2 in 4 charge bins (merged clusters)
Definition: SiPixelTemplate.h:162
SiPixelTemplateEntry::qmin
float qmin
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
Definition: SiPixelTemplate.h:123
SiPixelTemplateEntry::chi2xavgc2m
float chi2xavgc2m[4]
1st pass chi2 min search: average x chi^2 in 4 charge bins (merged clusters)
Definition: SiPixelTemplate.h:166
ysize
const Int_t ysize
Definition: trackSplitPlot.h:43
SiPixelTemplate::mpvvav2
float mpvvav2()
most probable charge in 2-cluster Vavilov distribution (not actually for larger kappa) ...
Definition: SiPixelTemplate.h:670
SiPixelTemplate::xflcorr
float xflcorr(int binq, float qflx)
Definition: SiPixelTemplate.cc:2448
SiPixelTemplateEntry::xgsiggen
float xgsiggen[4]
generic algorithm: average sigma_x from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:171
SiPixelTemplate::cotb_current_
float cotb_current_
current cot beta
Definition: SiPixelTemplate.h:692
cms::cuda::assert
assert(be >=bs)
SiPixelTemplateEntry::chi2xminc2m
float chi2xminc2m[4]
1st pass chi2 min search: minimum of x chi^2 in 4 charge bins (merged clusters)
Definition: SiPixelTemplate.h:167
SiPixelTemplateHeader::ID
int ID
&lt; template header structure
Definition: SiPixelTemplate.h:183
edm::isFinite
constexpr bool isFinite(T x)
SiPixelTemplateHeader::Bfield
float Bfield
Bfield in Tesla.
Definition: SiPixelTemplate.h:200
SiPixelTemplateEntry::cotalpha
float cotalpha
cot(alpha) is proportional to cluster length in x and is basis of interpolation
Definition: SiPixelTemplate.h:107
SiPixelTemplateHeader::temperature
float temperature
detector temperature in deg K
Definition: SiPixelTemplate.h:194
funct::arg
A arg
Definition: Factorize.h:31
SiPixelTemplateEntry::ytemp
float ytemp[9][21]
templates for y-reconstruction (binned over 1 central pixel)
Definition: SiPixelTemplate.h:139
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
SiPixelTemplate::sxone
float sxone()
rms for one pixel x-clusters
Definition: SiPixelTemplate.h:419
SiPixelTemplateEntry::yrmsc2m
float yrmsc2m[4]
1st pass chi2 min search: average y-rms of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:161
SiPixelTemplateHeader::lorywidth
float lorywidth
estimate of y-lorentz width for optimal resolution
Definition: SiPixelTemplate.h:189
BYM1
#define BYM1
Definition: SiPixelTemplateDefs.h:27
BXM2
#define BXM2
Definition: SiPixelTemplateDefs.h:36
SiPixelTemplateHeader::fbin
float fbin[3]
The QBin definitions in Q_clus/Q_avg.
Definition: SiPixelTemplate.h:201
SiPixelTemplate::ytemp
void ytemp(int fybin, int lybin, float ytemplate[41][21+4])
Definition: SiPixelTemplate.cc:2511
SiPixelTemplate::cxtemp
int cxtemp()
Return central pixel of x-template pixels above readout threshold.
Definition: SiPixelTemplate.cc:2709
SiPixelTemplateHeader::templ_version
int templ_version
Version number of the template to ensure code compatibility.
Definition: SiPixelTemplate.h:199
visualization-live-secondInstance_cfg.m
tuple m
Definition: visualization-live-secondInstance_cfg.py:79
SiPixelTemplateEntry::xpar
float xpar[2][5]
projected x-pixel uncertainty parameterization
Definition: SiPixelTemplate.h:140
edm::for
for(const auto &name:names)
Definition: PreMixingModule.cc:140
SiPixelTemplateEntry::ygx0
float ygx0[4]
average y0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:144
SiPixelTemplateEntry::xgx0gen
float xgx0gen[4]
generic algorithm: average x0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:170
SiPixelTemplate::kappavav2_
float kappavav2_
kappa parameter for 2-cluster Vavilov distribution
Definition: SiPixelTemplate.h:778
SiPixelTemplateEntry::symax
float symax
average pixel signal for y-projection of cluster
Definition: SiPixelTemplate.h:113
SiPixelTemplate::sigmavav_
float sigmavav_
&quot;sigma&quot; scale fctor for Vavilov distribution
Definition: SiPixelTemplate.h:774
SiPixelTemplateEntry::kappavav2
float kappavav2
kappa parameter for Vavilov distribution for 2 merged clusters
Definition: SiPixelTemplate.h:137
SiPixelTemplateEntry::xgsig
float xgsig[4]
average sigma_x from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:150
SiPixelTemplate::kappavav
float kappavav()
kappa parameter for Vavilov distribution
Definition: SiPixelTemplate.h:669
SiPixelTemplateEntry::cotbeta
float cotbeta
cot(beta) is proportional to cluster length in y and is basis of interpolation
Definition: SiPixelTemplate.h:109
SiPixelTemplateEntry::beta
float beta
beta track angle (defined in CMS CMS IN 2004/014)
Definition: SiPixelTemplate.h:108
SiPixelTemplate::vavilov_pars
void vavilov_pars(double &mpv, double &sigma, double &kappa)
Definition: SiPixelTemplate.cc:4002
DeadROC_duringRun.dir
string dir
Definition: DeadROC_duringRun.py:23
SiPixelTemplate::yrms
float yrms(int i)
average y-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:452
ENDL
#define ENDL
Definition: SiPixelTemplate.cc:113
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
SiPixelTemplate::enty0_
const SiPixelTemplateEntry * enty0_
Definition: SiPixelTemplate.h:799
SiPixelTemplate::lorywidth_
float lorywidth_
Lorentz y-width (sign corrected for fpix frame)
Definition: SiPixelTemplate.h:779
BHX
#define BHX
Definition: SiPixelTemplateDefs.h:34
if
if(conf_.getParameter< bool >("UseStripCablingDB"))
Definition: SiStripLAProfileBooker.cc:72
SiPixelTemplate::yflcorr
float yflcorr(int binq, float qfly)
Definition: SiPixelTemplate.cc:2392
SimplePixel::s
float s
distance from track entry
Definition: SimplePixel.h:24
SiPixelTemplateEntry::xrmsc2m
float xrmsc2m[4]
1st pass chi2 min search: average x-rms of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:165
change_name.diff
list diff
Definition: change_name.py:13
SiPixelTemplate::ytemp3d
void ytemp3d(int j, int k, std::vector< float > &ytemplate)
Definition: SiPixelTemplate.cc:2842
SiPixelTemplate::temperrors
void temperrors(int id, float cotalpha, float cotbeta, int qBin, float &sigmay, float &sigmax, float &sy1, float &sy2, float &sx1, float &sx2)
Definition: SiPixelTemplate.cc:3467
SiPixelTemplateHeader::lorxwidth
float lorxwidth
estimate of x-lorentz width for optimal resolution
Definition: SiPixelTemplate.h:190
SiPixelTemplateEntry::xgx0
float xgx0[4]
average x0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:149
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
SiPixelTemplate::ysigma2
void ysigma2(int fypix, int lypix, float sythr, float ysum[21+4], float ysig2[21+4])
SiPixelTemplateHeader::title
char title[80]
template title
Definition: SiPixelTemplate.h:198
SiPixelTemplate::sigmavav2_
float sigmavav2_
&quot;sigma&quot; scale fctor for 2-cluster Vavilov distribution
Definition: SiPixelTemplate.h:777
BYM2
#define BYM2
Definition: SiPixelTemplateDefs.h:28
SimplePixel::j
int j
y index of traversed pixel
Definition: SimplePixel.h:28
SiPixelTemplate::sigmavav2
float sigmavav2()
&quot;sigma&quot; scale fctor for 2-cluster Vavilov distribution
Definition: SiPixelTemplate.h:673
SiPixelTemplateEntry::ygsiggen
float ygsiggen[4]
generic algorithm: average sigma_y from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:169
SiPixelTemplateEntry::fracxone
float fracxone
fraction of sample with xsize = 1
Definition: SiPixelTemplate.h:174
SiPixelTemplate::sxtwo
float sxtwo()
rms for one double-pixel x-clusters
Definition: SiPixelTemplate.h:421
SiPixelTemplateEntry::mpvvav
float mpvvav
most probable charge in Vavilov distribution (not actually for larger kappa)
Definition: SiPixelTemplate.h:132
SiPixelTemplate::ysize_
float ysize_
Pixel y-size.
Definition: SiPixelTemplate.h:784
BYM3
#define BYM3
Definition: SiPixelTemplateDefs.h:29
SiPixelTemplate::xsize_
float xsize_
Pixel x-size.
Definition: SiPixelTemplate.h:783
SiPixelTemplate::kappavav_
float kappavav_
kappa parameter for Vavilov distribution
Definition: SiPixelTemplate.h:775
SiPixelTemplateEntry::sxtwo
float sxtwo
rms for one double-pixel x-clusters
Definition: SiPixelTemplate.h:122
SiPixelTemplateEntry::r_qMeas_qTrue
float r_qMeas_qTrue
ratio of measured to true cluster charge
Definition: SiPixelTemplate.h:178
SiPixelTemplateHeader::Dtype
int Dtype
detector type (0=BPix, 1=FPix)
Definition: SiPixelTemplate.h:187
SiPixelTemplateEntry::chi2yavg
float chi2yavg[4]
average y chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:152
SiPixelTemplateEntry::alpha
float alpha
alpha track angle (defined in CMS CMS IN 2004/014)
Definition: SiPixelTemplate.h:106
SiPixelTemplateEntry::xtemp
float xtemp[9][13]
templates for x-reconstruction (binned over 1 central pixel)
Definition: SiPixelTemplate.h:141
SiPixelTemplateStore::head
SiPixelTemplateHeader head
&lt; template storage structure
Definition: SiPixelTemplate.h:208
SiPixelTemplateEntry::chi2xminone
float chi2xminone
minimum of x chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:159
SiPixelTemplateEntry::ypar
float ypar[2][5]
projected y-pixel uncertainty parameterization
Definition: SiPixelTemplate.h:138
SiPixelTemplateStore::cotbetaX
std::array< float, 80 > cotbetaX
Definition: SiPixelTemplate.h:211
SiPixelTemplate::enty1_
const SiPixelTemplateEntry * enty1_
Definition: SiPixelTemplate.h:800
SiPixelTemplate::xtemp3d
void xtemp3d(int j, int k, std::vector< float > &xtemplate)
Definition: SiPixelTemplate.cc:2944
SiPixelTemplate::thePixelTemp_
const std::vector< SiPixelTemplateStore > & thePixelTemp_
Definition: SiPixelTemplate.h:812
SiPixelTemplateHeader::s50
float s50
1/2 of the multihit dcol threshold in electrons
Definition: SiPixelTemplate.h:196
TYSIZE
#define TYSIZE
Definition: SiPixelTemplateDefs.h:21
SiPixelTemplate::xtemp
void xtemp(int fxbin, int lxbin, float xtemplate[41][13+4])
Definition: SiPixelTemplate.cc:2587
SiPixelTemplateHeader::zsize
float zsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate.h:204
SiPixelTemplate::qbin_dist
void qbin_dist(int id, float cotalpha, float cotbeta, float qbin_frac[4], float &ny1_frac, float &ny2_frac, float &nx1_frac, float &nx2_frac)
Definition: SiPixelTemplate.cc:3627
SiPixelTemplateHeader::Vbias
float Vbias
detector bias potential in Volts
Definition: SiPixelTemplate.h:193
SiPixelTemplateEntry::qbfrac
float qbfrac[3]
fraction of sample in qbin = 0-2 (&gt;=3 is the complement)
Definition: SiPixelTemplate.h:172
SiPixelTemplateHeader::xsize
float xsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate.h:202
SiPixelTemplate::ytemp3d_int
void ytemp3d_int(int nypix, int &nybins)
Definition: SiPixelTemplate.cc:2764
SiPixelTemplate::qavg_avg_
float qavg_avg_
average of cluster charge less than qavg
Definition: SiPixelTemplate.h:786
SiPixelTemplateEntry::fracyone
float fracyone
fraction of sample with ysize = 1
Definition: SiPixelTemplate.h:173
SiPixelTemplateEntry::qmin2
float qmin2
tighter minimum cluster charge for valid hit (keeps 99.8% of simulated hits)
Definition: SiPixelTemplate.h:124
SiPixelTemplate::yxratio
float yxratio()
fractional distance in y between cotalpha templates slices
Definition: SiPixelTemplate.h:440
SiPixelTemplateEntry::yavg
float yavg[4]
average y-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:142
SiPixelTemplate::xtemp3d_int
void xtemp3d_int(int nxpix, int &nxbins)
Definition: SiPixelTemplate.cc:2865
BHY
#define BHY
Definition: SiPixelTemplateDefs.h:26
SiPixelTemplate::pushfile
static bool pushfile(int filenum, std::vector< SiPixelTemplateStore > &pixelTemp, std::string dir="CalibTracker/SiPixelESProducers/data/")
Definition: SiPixelTemplate.cc:130
SiPixelTemplateEntry::chi2ymin
float chi2ymin[4]
minimum of y chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:153
SiPixelTemplateEntry::yavgc2m
float yavgc2m[4]
1st pass chi2 min search: average y-bias of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:160
SimplePixel::x
float x
x coordinate of boundary intersection
Definition: SimplePixel.h:25
SiPixelTemplateHeader::NTy
int NTy
number of Template y entries
Definition: SiPixelTemplate.h:184
SiPixelTemplateEntry::dxtwo
float dxtwo
mean offset/correction for one double-pixel x-clusters
Definition: SiPixelTemplate.h:121
SiPixelTemplateHeader::lorxbias
float lorxbias
estimate of x-lorentz bias
Definition: SiPixelTemplate.h:192
cuda_std::lower_bound
__host__ __device__ constexpr RandomIt lower_bound(RandomIt first, RandomIt last, const T &value, Compare comp={})
Definition: cudastdAlgorithm.h:27
SiPixelTemplateStore::entx
SiPixelTemplateEntry entx[80][80]
Definition: SiPixelTemplate.h:216
LOGINFO
#define LOGINFO(x)
Definition: SiPixelTemplate.cc:111
SimplePixel::i
int i
x index of traversed pixel
Definition: SimplePixel.h:27
SiPixelTemplate::sigmavav
float sigmavav()
&quot;sigma&quot; scale fctor for Vavilov distribution
Definition: SiPixelTemplate.h:668
SiPixelTemplateHeader::NTyx
int NTyx
number of Template y-slices of x entries
Definition: SiPixelTemplate.h:185
mps_splice.entry
list entry
Definition: mps_splice.py:68
SiPixelTemplateEntry::yflpar
float yflpar[4][6]
Aqfl-parameterized y-correction in 4 charge bins.
Definition: SiPixelTemplate.h:146
edm::FileInPath::fullPath
std::string fullPath() const
Definition: FileInPath.cc:161
SiPixelTemplateEntry::sxone
float sxone
rms for one pixel x-clusters
Definition: SiPixelTemplate.h:118
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144
SiPixelTemplateEntry::chi2xmin
float chi2xmin[4]
minimum of x chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:155
SiPixelTemplate::lorxwidth_
float lorxwidth_
Lorentz x-width.
Definition: SiPixelTemplate.h:780
SiPixelTemplate::index_id_
int index_id_
current index
Definition: SiPixelTemplate.h:690
SiPixelTemplateEntry::chi2yminone
float chi2yminone
minimum of y chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:157
SiPixelTemplateStore::enty
SiPixelTemplateEntry enty[100]
60 Barrel x templates spanning cluster lengths from -6px (-1.125Rad) to +6px (+1.125Rad) in each of 6...
Definition: SiPixelTemplate.h:214
SiPixelTemplate::qbin
int qbin(int id, float cotalpha, float cotbeta, float locBz, float locBx, float qclus, float &pixmx, float &sigmay, float &deltay, float &sigmax, float &deltax, float &sy1, float &dy1, float &sy2, float &dy2, float &sx1, float &dx1, float &sx2, float &dx2)
Definition: SiPixelTemplate.cc:2991
PixelTestBeamValidation_cfi.Ny
int Ny
Definition: PixelTestBeamValidation_cfi.py:8
SiPixelTemplateEntry::costrk
float costrk[3]
direction cosines of tracks used to generate this entry
Definition: SiPixelTemplate.h:110
SiPixelTemplate::postInit
static void postInit(std::vector< SiPixelTemplateStore > &thePixelTemp_)
Definition: SiPixelTemplate.cc:1288
SiPixelTemplate::mpvvav2_
float mpvvav2_
most probable charge in 2-cluster Vavilov distribution (not actually for larger kappa) ...
Definition: SiPixelTemplate.h:776
SiPixelTemplateEntry::ygsig
float ygsig[4]
average sigma_y from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:145
SiPixelTemplateEntry::sigmavav2
float sigmavav2
&quot;sigma&quot; scale fctor for Vavilov distribution for 2 merged clusters
Definition: SiPixelTemplate.h:136
SiPixelTemplateEntry::sytwo
float sytwo
rms for one double-pixel y-clusters
Definition: SiPixelTemplate.h:120
xsize
const Int_t xsize
Definition: trackSplitPlot.h:42
SiPixelTemplateEntry::xrmsgen
float xrmsgen[4]
generic algorithm: average x-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:128
SiPixelTemplate::id_current_
int id_current_
current id
Definition: SiPixelTemplate.h:689
SiPixelTemplate::zsize_
float zsize_
Pixel z-size (thickness)
Definition: SiPixelTemplate.h:785
SiPixelTemplateEntry::chi2xavg
float chi2xavg[4]
average x chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:154
SiPixelTemplate::yratio_
float yratio_
fractional distance in y between cotbeta templates
Definition: SiPixelTemplate.h:730
SiPixelTemplate::xxratio
float xxratio()
fractional distance in x between cotalpha templates
Definition: SiPixelTemplate.h:441
SiPixelTemplateEntry::kappavav
float kappavav
kappa parameter for Vavilov distribution
Definition: SiPixelTemplate.h:134
SiPixelTemplateEntry::xavggen
float xavggen[4]
generic algorithm: average x-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:127
SiPixelTemplateEntry::yrmsgen
float yrmsgen[4]
generic algorithm: average y-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:126
SiPixelTemplateEntry::clsleny
float clsleny
cluster y-length in pixels at signal height symax/2
Definition: SiPixelTemplate.h:130
Generated for CMSSW Reference Manual by doxygen 1.8.5