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