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