CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
SiPixelTemplate2D.cc
Go to the documentation of this file.
1 //
2 // SiPixelTemplate2D.cc Version 2.65
3 //
4 // Full 2-D templates for cluster splitting, simulated cluster reweighting, and improved cluster probability
5 //
6 // Created by Morris Swartz on 12/01/09.
7 // 2009 __TheJohnsHopkinsUniversity__.
8 //
9 // V1.01 - fix qavg_ filling
10 // V1.02 - Add locBz to test if FPix use is out of range
11 // V1.03 - Fix edge checking on final template to increase template size and to properly truncate cluster
12 // v2.00 - Major changes to accommodate 2D reconstruction
13 // v2.10 - Change chi2 and error scaling information to work with partially reconstructed clusters
14 // v2.20 - Add cluster charge probability information, side loading for template generation
15 // v2.21 - Double derivative interval [improves fit convergence]
16 // v2.25 - Resize template store to accommodate FPix Templates
17 // v2.30 - Fix bug found by P. Shuetze that compromises sqlite file loading
18 // v2.35 - Add directory path selection to the ascii pushfile method
19 // v2.50 - Change template storage to dynamically allocated 2D arrays of SiPixelTemplateEntry2D structs
20 // v2.51 - Ensure that the derivative arrays are correctly zeroed between calls
21 // v2.52 - Improve cosmetics for increased style points from judges
22 // v2.60 - Fix FPix multiframe lookup problem [takes +-cotalpha and +-cotbeta]
23 // v2.61a - Code 2.60 fix correctly
24 // v2.65 - change double pixel flags to work with new shifted reco code definition
25 //
26 
27 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
28 #else
29 #include <math.h>
30 #endif
31 #include <algorithm>
32 #include <vector>
33 //#include "boost/multi_array.hpp"
34 #include <iostream>
35 #include <iomanip>
36 #include <sstream>
37 #include <fstream>
38 
39 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
43 #define LOGERROR(x) LogError(x)
44 #define LOGINFO(x) LogInfo(x)
45 #define ENDL " "
47 using namespace edm;
48 #else
49 #include "SiPixelTemplate2D.h"
50 #define LOGERROR(x) std::cout << x << ": "
51 #define LOGINFO(x) std::cout << x << ": "
52 #define ENDL std::endl
53 #endif
54 
55 //****************************************************************
60 //****************************************************************
61 bool SiPixelTemplate2D::pushfile(int filenum, std::vector<SiPixelTemplateStore2D>& pixelTemp, std::string dir) {
62  // Add template stored in external file numbered filenum to theTemplateStore
63 
64  // Local variables
65  const int code_version = {21};
66 
67  // Create a filename for this run
68  std::string tempfile = std::to_string(filenum);
69 
70  // Create different path in CMSSW than standalone
71 
72 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
73  // If integer filenum has less than 4 digits, prepend 0's until we have four numerical characters, e.g. "0292"
74  int nzeros = 4 - tempfile.length();
75  if (nzeros > 0)
76  tempfile = std::string(nzeros, '0') + tempfile;
78 
79  tempfile = dir + "template_summary2D_zp" + tempfile + ".out";
80  edm::FileInPath file(tempfile); // Find the file in CMSSW
81  tempfile = file.fullPath(); // Put it back with the whole path.
82 
83 #else
84  // This is the same as above, but more elegant.
85  std::ostringstream tout;
86  tout << "template_summary2D_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
87  tempfile = tout.str();
88 
89 #endif
90 
91  // Open the template file
92  //
93  std::ifstream in_file(tempfile);
94  if (in_file.is_open() && in_file.good()) {
95  // Create a local template storage entry
96  SiPixelTemplateStore2D theCurrentTemp;
97 
98  // Read-in a header string first and print it
99  char c;
100  for (int i = 0; (c = in_file.get()) != '\n'; ++i) {
101  if (i < 79) {
102  theCurrentTemp.head.title[i] = c;
103  } else {
104  theCurrentTemp.head.title[79] = '\0';
105  }
106  }
107  LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
108 
109  // next, the header information
110  in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
111  theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
112  theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
113  theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
114  theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
115  theCurrentTemp.head.zsize;
116 
117  if (in_file.fail()) {
118  LOGERROR("SiPixelTemplate2D") << "Error reading file 0A, no template load" << ENDL;
119  return false;
120  }
121 
122  if (theCurrentTemp.head.templ_version > 17) {
123  in_file >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
124  theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
125 
126  if (in_file.fail()) {
127  LOGERROR("SiPixelTemplate2D") << "Error reading file 0B, no template load" << ENDL;
128  return false;
129  }
130  } else {
131  // This is for older [legacy] payloads
132  theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
133  theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth / 2.f;
134  theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth / 2.f;
135  theCurrentTemp.head.fbin[0] = 1.5f;
136  theCurrentTemp.head.fbin[1] = 1.00f;
137  theCurrentTemp.head.fbin[2] = 0.85f;
138  }
139 
140  LOGINFO("SiPixelTemplate2D") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
141  << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
142  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
143  << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
144  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
145  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
146  << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
147  << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold "
148  << theCurrentTemp.head.ss50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth
149  << ", y Lorentz Bias " << theCurrentTemp.head.lorybias << ", x Lorentz width "
150  << theCurrentTemp.head.lorxwidth << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
151  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
152  << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
153  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
154  << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
155 
156  if (theCurrentTemp.head.templ_version < code_version) {
157  LOGERROR("SiPixelTemplate2D") << "code expects version " << code_version << ", no template load" << ENDL;
158  return false;
159  }
160 
161  if (theCurrentTemp.head.NTy != 0) {
162  LOGERROR("SiPixelTemplate2D")
163  << "Trying to load 1-d template info into the 2-d template object, check your DB/global tag!" << ENDL;
164  return false;
165  }
166 
167  // next, layout the 2-d structure needed to store template
168 
169  theCurrentTemp.entry.resize(theCurrentTemp.head.NTyx);
170  for (auto& item : theCurrentTemp.entry)
171  item.resize(theCurrentTemp.head.NTxx);
172 
173  // Read in the file info
174 
175  for (int iy = 0; iy < theCurrentTemp.head.NTyx; ++iy) {
176  for (int jx = 0; jx < theCurrentTemp.head.NTxx; ++jx) {
177  in_file >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0] >>
178  theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2];
179 
180  if (in_file.fail()) {
181  LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # "
182  << theCurrentTemp.entry[iy][jx].runnum << ENDL;
183  return false;
184  }
185 
186  // Calculate cot(alpha) and cot(beta) for this entry
187 
188  theCurrentTemp.entry[iy][jx].cotalpha =
189  theCurrentTemp.entry[iy][jx].costrk[0] / theCurrentTemp.entry[iy][jx].costrk[2];
190 
191  theCurrentTemp.entry[iy][jx].cotbeta =
192  theCurrentTemp.entry[iy][jx].costrk[1] / theCurrentTemp.entry[iy][jx].costrk[2];
193 
194  in_file >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >>
195  theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin >>
196  theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >>
197  theCurrentTemp.entry[iy][jx].jxmax;
198 
199  if (in_file.fail()) {
200  LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # "
201  << theCurrentTemp.entry[iy][jx].runnum << ENDL;
202  return false;
203  }
204 
205  for (int k = 0; k < 2; ++k) {
206  in_file >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1] >>
207  theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >>
208  theCurrentTemp.entry[iy][jx].xypar[k][4];
209 
210  if (in_file.fail()) {
211  LOGERROR("SiPixelTemplate2D")
212  << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
213  return false;
214  }
215  }
216 
217  for (int k = 0; k < 2; ++k) {
218  in_file >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1] >>
219  theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >>
220  theCurrentTemp.entry[iy][jx].lanpar[k][4];
221 
222  if (in_file.fail()) {
223  LOGERROR("SiPixelTemplate2D")
224  << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
225  return false;
226  }
227  }
228 
229  // Read the 2D template entries as floats [they are formatted that way] and cast to short ints
230 
231  float dummy[T2YSIZE];
232  for (int l = 0; l < 7; ++l) {
233  for (int k = 0; k < 7; ++k) {
234  for (int j = 0; j < T2XSIZE; ++j) {
235  for (int i = 0; i < T2YSIZE; ++i) {
236  in_file >> dummy[i];
237  }
238  if (in_file.fail()) {
239  LOGERROR("SiPixelTemplate2D")
240  << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
241  return false;
242  }
243  for (int i = 0; i < T2YSIZE; ++i) {
244  theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j] = (short int)dummy[i];
245  }
246  }
247  }
248  }
249 
250  in_file >> theCurrentTemp.entry[iy][jx].chi2ppix >> theCurrentTemp.entry[iy][jx].chi2scale >>
251  theCurrentTemp.entry[iy][jx].offsetx[0] >> theCurrentTemp.entry[iy][jx].offsetx[1] >>
252  theCurrentTemp.entry[iy][jx].offsetx[2] >> theCurrentTemp.entry[iy][jx].offsetx[3] >>
253  theCurrentTemp.entry[iy][jx].offsety[0] >> theCurrentTemp.entry[iy][jx].offsety[1] >>
254  theCurrentTemp.entry[iy][jx].offsety[2] >> theCurrentTemp.entry[iy][jx].offsety[3];
255 
256  if (in_file.fail()) {
257  LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # "
258  << theCurrentTemp.entry[iy][jx].runnum << ENDL;
259  return false;
260  }
261 
262  in_file >> theCurrentTemp.entry[iy][jx].clsleny >> theCurrentTemp.entry[iy][jx].clslenx >>
263  theCurrentTemp.entry[iy][jx].mpvvav >> theCurrentTemp.entry[iy][jx].sigmavav >>
264  theCurrentTemp.entry[iy][jx].kappavav >> theCurrentTemp.entry[iy][jx].scalexavg >>
265  theCurrentTemp.entry[iy][jx].scaleyavg >> theCurrentTemp.entry[iy][jx].delyavg >>
266  theCurrentTemp.entry[iy][jx].delysig >> theCurrentTemp.entry[iy][jx].spare[0];
267 
268  if (in_file.fail()) {
269  LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # "
270  << theCurrentTemp.entry[iy][jx].runnum << ENDL;
271  return false;
272  }
273 
274  in_file >> theCurrentTemp.entry[iy][jx].scalex[0] >> theCurrentTemp.entry[iy][jx].scalex[1] >>
275  theCurrentTemp.entry[iy][jx].scalex[2] >> theCurrentTemp.entry[iy][jx].scalex[3] >>
276  theCurrentTemp.entry[iy][jx].scaley[0] >> theCurrentTemp.entry[iy][jx].scaley[1] >>
277  theCurrentTemp.entry[iy][jx].scaley[2] >> theCurrentTemp.entry[iy][jx].scaley[3] >>
278  theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2];
279 
280  if (in_file.fail()) {
281  LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # "
282  << theCurrentTemp.entry[iy][jx].runnum << ENDL;
283  return false;
284  }
285  }
286  }
287 
288  in_file.close();
289 
290  // Add this template to the store
291 
292  pixelTemp.push_back(theCurrentTemp);
293 
294  return true;
295 
296  } else {
297  // If file didn't open, report this
298 
299  LOGERROR("SiPixelTemplate2D") << "Error opening File" << tempfile << ENDL;
300  return false;
301  }
302 
303 } // TempInit
304 
305 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
306 
307 //****************************************************************
311 //****************************************************************
313  std::vector<SiPixelTemplateStore2D>& pixelTemp) {
314  // Add template stored in external dbobject to theTemplateStore
315 
316  const int code_version = {21};
317 
318  // We must create a new object because dbobject must be a const and our stream must not be
319  SiPixel2DTemplateDBObject db = dbobject;
320 
321  // Create a local template storage entry
322  SiPixelTemplateStore2D theCurrentTemp;
323 
324  // Fill the template storage for each template calibration stored in the db
325  for (int m = 0; m < db.numOfTempl(); ++m) {
326  // Read-in a header string first and print it
327 
329  for (int i = 0; i < 20; ++i) {
330  temp.f = db.sVector()[db.index()];
331  theCurrentTemp.head.title[4 * i] = temp.c[0];
332  theCurrentTemp.head.title[4 * i + 1] = temp.c[1];
333  theCurrentTemp.head.title[4 * i + 2] = temp.c[2];
334  theCurrentTemp.head.title[4 * i + 3] = temp.c[3];
335  db.incrementIndex(1);
336  }
337  theCurrentTemp.head.title[79] = '\0';
338  LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
339 
340  // next, the header information
341 
342  db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >>
343  theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >>
344  theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >>
345  theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >>
346  theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
347  theCurrentTemp.head.zsize;
348 
349  if (db.fail()) {
350  LOGERROR("SiPixelTemplate2D") << "Error reading file 0A, no template load" << ENDL;
351  return false;
352  }
353 
354  LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title
355  << " code version = " << code_version << " object version "
356  << theCurrentTemp.head.templ_version << ENDL;
357 
358  if (theCurrentTemp.head.templ_version > 17) {
359  db >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
360  theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
361 
362  if (db.fail()) {
363  LOGERROR("SiPixelTemplate2D") << "Error reading file 0B, no template load" << ENDL;
364  return false;
365  }
366  } else {
367  // This is for older [legacy] payloads and the numbers are indeed magic [they are part of the payload for v>17]
368  theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
369  theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth / 2.f;
370  theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth / 2.f;
371  theCurrentTemp.head.fbin[0] = 1.50f;
372  theCurrentTemp.head.fbin[1] = 1.00f;
373  theCurrentTemp.head.fbin[2] = 0.85f;
374  }
375 
376  LOGINFO("SiPixelTemplate2D") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
377  << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
378  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx
379  << ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
380  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
381  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
382  << ", Q-scaling factor " << theCurrentTemp.head.qscale << ", 1/2 multi dcol threshold "
383  << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold "
384  << theCurrentTemp.head.ss50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth
385  << ", y Lorentz Bias " << theCurrentTemp.head.lorybias << ", x Lorentz width "
386  << theCurrentTemp.head.lorxwidth << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
387  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", "
388  << theCurrentTemp.head.fbin[1] << ", " << theCurrentTemp.head.fbin[2]
389  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size "
390  << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
391 
392  if (theCurrentTemp.head.templ_version < code_version) {
393  LOGINFO("SiPixelTemplate2D") << "code expects version " << code_version << " finds "
394  << theCurrentTemp.head.templ_version << ", load anyway " << ENDL;
395  }
396 
397  if (theCurrentTemp.head.NTy != 0) {
398  LOGERROR("SiPixelTemplate2D")
399  << "Trying to load 1-d template info into the 2-d template object, check your DB/global tag!" << ENDL;
400  return false;
401  }
402 
403  // next, layout the 2-d structure needed to store template
404 
405  theCurrentTemp.entry.resize(theCurrentTemp.head.NTyx);
406  for (auto& item : theCurrentTemp.entry)
407  item.resize(theCurrentTemp.head.NTxx);
408 
409  // Read in the file info
410 
411  for (int iy = 0; iy < theCurrentTemp.head.NTyx; ++iy) {
412  for (int jx = 0; jx < theCurrentTemp.head.NTxx; ++jx) {
413  db >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0] >>
414  theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2];
415 
416  if (db.fail()) {
417  LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # "
418  << theCurrentTemp.entry[iy][jx].runnum << ENDL;
419  return false;
420  }
421 
422  // Calculate cot(alpha) and cot(beta) for this entry
423 
424  theCurrentTemp.entry[iy][jx].cotalpha =
425  theCurrentTemp.entry[iy][jx].costrk[0] / theCurrentTemp.entry[iy][jx].costrk[2];
426 
427  theCurrentTemp.entry[iy][jx].cotbeta =
428  theCurrentTemp.entry[iy][jx].costrk[1] / theCurrentTemp.entry[iy][jx].costrk[2];
429 
430  db >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >>
431  theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin >>
432  theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >>
433  theCurrentTemp.entry[iy][jx].jxmax;
434 
435  if (db.fail()) {
436  LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # "
437  << theCurrentTemp.entry[iy][jx].runnum << ENDL;
438  return false;
439  }
440 
441  for (int k = 0; k < 2; ++k) {
442  db >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1] >>
443  theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >>
444  theCurrentTemp.entry[iy][jx].xypar[k][4];
445 
446  if (db.fail()) {
447  LOGERROR("SiPixelTemplate2D")
448  << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
449  return false;
450  }
451  }
452 
453  for (int k = 0; k < 2; ++k) {
454  db >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1] >>
455  theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >>
456  theCurrentTemp.entry[iy][jx].lanpar[k][4];
457 
458  if (db.fail()) {
459  LOGERROR("SiPixelTemplate2D")
460  << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
461  return false;
462  }
463  }
464 
465  // Read the 2D template entries as floats [they are formatted that way] and cast to short ints
466 
467  float dummy[T2YSIZE];
468  for (int l = 0; l < 7; ++l) {
469  for (int k = 0; k < 7; ++k) {
470  for (int j = 0; j < T2XSIZE; ++j) {
471  for (int i = 0; i < T2YSIZE; ++i) {
472  db >> dummy[i];
473  }
474  if (db.fail()) {
475  LOGERROR("SiPixelTemplate2D")
476  << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
477  return false;
478  }
479  for (int i = 0; i < T2YSIZE; ++i) {
480  theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j] = (short int)dummy[i];
481  }
482  }
483  }
484  }
485 
486  db >> theCurrentTemp.entry[iy][jx].chi2ppix >> theCurrentTemp.entry[iy][jx].chi2scale >>
487  theCurrentTemp.entry[iy][jx].offsetx[0] >> theCurrentTemp.entry[iy][jx].offsetx[1] >>
488  theCurrentTemp.entry[iy][jx].offsetx[2] >> theCurrentTemp.entry[iy][jx].offsetx[3] >>
489  theCurrentTemp.entry[iy][jx].offsety[0] >> theCurrentTemp.entry[iy][jx].offsety[1] >>
490  theCurrentTemp.entry[iy][jx].offsety[2] >> theCurrentTemp.entry[iy][jx].offsety[3];
491 
492  if (db.fail()) {
493  LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # "
494  << theCurrentTemp.entry[iy][jx].runnum << ENDL;
495  return false;
496  }
497 
498  db >> theCurrentTemp.entry[iy][jx].clsleny >> theCurrentTemp.entry[iy][jx].clslenx >>
499  theCurrentTemp.entry[iy][jx].mpvvav >> theCurrentTemp.entry[iy][jx].sigmavav >>
500  theCurrentTemp.entry[iy][jx].kappavav >> theCurrentTemp.entry[iy][jx].scalexavg >>
501  theCurrentTemp.entry[iy][jx].scaleyavg >> theCurrentTemp.entry[iy][jx].delyavg >>
502  theCurrentTemp.entry[iy][jx].delysig >> theCurrentTemp.entry[iy][jx].spare[0];
503 
504  if (db.fail()) {
505  LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # "
506  << theCurrentTemp.entry[iy][jx].runnum << ENDL;
507  return false;
508  }
509 
510  db >> theCurrentTemp.entry[iy][jx].scalex[0] >> theCurrentTemp.entry[iy][jx].scalex[1] >>
511  theCurrentTemp.entry[iy][jx].scalex[2] >> theCurrentTemp.entry[iy][jx].scalex[3] >>
512  theCurrentTemp.entry[iy][jx].scaley[0] >> theCurrentTemp.entry[iy][jx].scaley[1] >>
513  theCurrentTemp.entry[iy][jx].scaley[2] >> theCurrentTemp.entry[iy][jx].scaley[3] >>
514  theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2];
515 
516  if (db.fail()) {
517  LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # "
518  << theCurrentTemp.entry[iy][jx].runnum << ENDL;
519  return false;
520  }
521  }
522  }
523 
524  // Add this template to the store
525 
526  pixelTemp.push_back(theCurrentTemp);
527  }
528 
529  return true;
530 
531 } // TempInit
532 
533 #endif
534 
536  if (id != id_current_) {
537  // Find the index corresponding to id
538 
539  index_id_ = -1;
540  for (int i = 0; i < (int)thePixelTemp_.size(); ++i) {
541  if (id == thePixelTemp_[i].head.ID) {
542  index_id_ = i;
543  id_current_ = id;
544 
545  // Copy the detector type to the private variable
546 
547  Dtype_ = thePixelTemp_[index_id_].head.Dtype;
548 
549  // Copy the charge scaling factor to the private variable
550 
551  qscale_ = thePixelTemp_[index_id_].head.qscale;
552 
553  // Copy the pseudopixel signal size to the private variable
554 
555  s50_ = thePixelTemp_[index_id_].head.s50;
556 
557  // Copy Qbinning info to private variables
558 
559  for (int j = 0; j < 3; ++j) {
560  fbin_[j] = thePixelTemp_[index_id_].head.fbin[j];
561  }
562 
563  // Copy the Lorentz widths to private variables
564 
565  lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
566  lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
567 
568  // Copy the pixel sizes private variables
569 
570  xsize_ = thePixelTemp_[index_id_].head.xsize;
571  ysize_ = thePixelTemp_[index_id_].head.ysize;
572  zsize_ = thePixelTemp_[index_id_].head.zsize;
573 
574  // Determine the size of this template
575 
576  Nyx_ = thePixelTemp_[index_id_].head.NTyx;
577  Nxx_ = thePixelTemp_[index_id_].head.NTxx;
578 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
579  if (Nyx_ < 2 || Nxx_ < 2) {
580  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_
581  << "has too few entries: Nyx/Nxx = " << Nyx_ << "/" << Nxx_ << std::endl;
582  }
583 #else
584  assert(Nyx_ > 1 && Nxx_ > 1);
585 #endif
586  int imidx = Nxx_ / 2;
587 
588  cotalpha0_ = thePixelTemp_[index_id_].entry[0][0].cotalpha;
589  cotalpha1_ = thePixelTemp_[index_id_].entry[0][Nxx_ - 1].cotalpha;
590  deltacota_ = (cotalpha1_ - cotalpha0_) / (float)(Nxx_ - 1);
591 
592  cotbeta0_ = thePixelTemp_[index_id_].entry[0][imidx].cotbeta;
593  cotbeta1_ = thePixelTemp_[index_id_].entry[Nyx_ - 1][imidx].cotbeta;
594  deltacotb_ = (cotbeta1_ - cotbeta0_) / (float)(Nyx_ - 1);
595 
596  break;
597  }
598  }
599  }
600 
601 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
602  if (index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
603  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
604  << ", Are you using the correct global tag?" << std::endl;
605  }
606 #else
607  assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
608 #endif
609  return true;
610 }
611 
612 // *************************************************************************************************************************************
623 // *************************************************************************************************************************************
624 
625 bool SiPixelTemplate2D::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx) {
626  // Interpolate for a new set of track angles
627 
628  // Local variables
629 
630  float acotb, dcota, dcotb;
631 
632  // Check to see if interpolation is valid
633 
634  if (id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
635  cota_current_ = cotalpha;
636  cotb_current_ = cotbeta;
637  // Try to find the correct template. Fill the class variable index_id_ .
638  success_ = getid(id);
639  }
640 
641 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
642  if (index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
643  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
644  << ", Are you using the correct global tag?" << std::endl;
645  }
646 #else
647  assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
648 #endif
649 
650  // Check angle limits and et up interpolation parameters
651 
652  float cota = cotalpha;
653  flip_x_ = false;
654  flip_y_ = false;
655  switch (Dtype_) {
656  case 0:
657  if (cotbeta < 0.f) {
658  flip_y_ = true;
659  }
660  break;
661  case 1:
662  if (locBz > 0.f) {
663  flip_y_ = true;
664  }
665  break;
666  case 2:
667  case 3:
668  case 4:
669  case 5:
670  if (locBx * locBz < 0.f) {
671  cota = std::abs(cotalpha);
672  flip_x_ = true;
673  }
674  if (locBx < 0.f) {
675  flip_y_ = true;
676  }
677  break;
678  default:
679 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
680  throw cms::Exception("DataCorrupt")
681  << "SiPixelTemplate2D::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
682 #else
683  std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
684 #endif
685  }
686 
687  if (cota < cotalpha0_) {
688  success_ = false;
689  jx0_ = 0;
690  jx1_ = 1;
691  adcota_ = 0.f;
692  } else if (cota > cotalpha1_) {
693  success_ = false;
694  jx0_ = Nxx_ - 1;
695  jx1_ = jx0_ - 1;
696  adcota_ = 0.f;
697  } else {
698  jx0_ = (int)((cota - cotalpha0_) / deltacota_ + 0.5f);
699  dcota = (cota - (cotalpha0_ + jx0_ * deltacota_)) / deltacota_;
700  adcota_ = fabs(dcota);
701  if (dcota > 0.f) {
702  jx1_ = jx0_ + 1;
703  if (jx1_ > Nxx_ - 1)
704  jx1_ = jx0_ - 1;
705  } else {
706  jx1_ = jx0_ - 1;
707  if (jx1_ < 0)
708  jx1_ = jx0_ + 1;
709  }
710  }
711 
712  // Interpolate the absolute value of cot(beta)
713 
714  acotb = std::abs(cotbeta);
715 
716  if (acotb < cotbeta0_) {
717  success_ = false;
718  iy0_ = 0;
719  iy1_ = 1;
720  adcotb_ = 0.f;
721  } else if (acotb > cotbeta1_) {
722  success_ = false;
723  iy0_ = Nyx_ - 1;
724  iy1_ = iy0_ - 1;
725  adcotb_ = 0.f;
726  } else {
727  iy0_ = (int)((acotb - cotbeta0_) / deltacotb_ + 0.5f);
728  dcotb = (acotb - (cotbeta0_ + iy0_ * deltacotb_)) / deltacotb_;
729  adcotb_ = fabs(dcotb);
730  if (dcotb > 0.f) {
731  iy1_ = iy0_ + 1;
732  if (iy1_ > Nyx_ - 1)
733  iy1_ = iy0_ - 1;
734  } else {
735  iy1_ = iy0_ - 1;
736  if (iy1_ < 0)
737  iy1_ = iy0_ + 1;
738  }
739  }
740 
741  // Calculate signed quantities
742 
743  lorydrift_ = lorywidth_ / 2.;
744  if (flip_y_)
745  lorydrift_ = -lorydrift_;
746  lorxdrift_ = lorxwidth_ / 2.;
747  if (flip_x_)
748  lorxdrift_ = -lorxdrift_;
749 
750  // Use pointers to the three angle pairs used in the interpolation
751 
752  entry00_ = &thePixelTemp_[index_id_].entry[iy0_][jx0_];
753  entry10_ = &thePixelTemp_[index_id_].entry[iy1_][jx0_];
754  entry01_ = &thePixelTemp_[index_id_].entry[iy0_][jx1_];
755 
756  // Interpolate things in cot(alpha)-cot(beta)
757 
758  qavg_ = entry00_->qavg + adcota_ * (entry01_->qavg - entry00_->qavg) + adcotb_ * (entry10_->qavg - entry00_->qavg);
759 
760  pixmax_ = entry00_->pixmax + adcota_ * (entry01_->pixmax - entry00_->pixmax) +
761  adcotb_ * (entry10_->pixmax - entry00_->pixmax);
762 
763  sxymax_ = entry00_->sxymax + adcota_ * (entry01_->sxymax - entry00_->sxymax) +
764  adcotb_ * (entry10_->sxymax - entry00_->sxymax);
765 
766  chi2avgone_ = entry00_->chi2avgone + adcota_ * (entry01_->chi2avgone - entry00_->chi2avgone) +
767  adcotb_ * (entry10_->chi2avgone - entry00_->chi2avgone);
768 
769  chi2minone_ = entry00_->chi2minone + adcota_ * (entry01_->chi2minone - entry00_->chi2minone) +
770  adcotb_ * (entry10_->chi2minone - entry00_->chi2minone);
771 
772  clsleny_ = entry00_->clsleny + adcota_ * (entry01_->clsleny - entry00_->clsleny) +
773  adcotb_ * (entry10_->clsleny - entry00_->clsleny);
774 
775  clslenx_ = entry00_->clslenx + adcota_ * (entry01_->clslenx - entry00_->clslenx) +
776  adcotb_ * (entry10_->clslenx - entry00_->clslenx);
777 
778  chi2ppix_ = entry00_->chi2ppix + adcota_ * (entry01_->chi2ppix - entry00_->chi2ppix) +
779  adcotb_ * (entry10_->chi2ppix - entry00_->chi2ppix);
780 
781  chi2scale_ = entry00_->chi2scale + adcota_ * (entry01_->chi2scale - entry00_->chi2scale) +
782  adcotb_ * (entry10_->chi2scale - entry00_->chi2scale);
783 
784  scaleyavg_ = entry00_->scaleyavg + adcota_ * (entry01_->scaleyavg - entry00_->scaleyavg) +
785  adcotb_ * (entry10_->scaleyavg - entry00_->scaleyavg);
786 
787  scalexavg_ = entry00_->scalexavg + adcota_ * (entry01_->scalexavg - entry00_->scalexavg) +
788  adcotb_ * (entry10_->scalexavg - entry00_->scalexavg);
789 
790  delyavg_ = entry00_->delyavg + adcota_ * (entry01_->delyavg - entry00_->delyavg) +
791  adcotb_ * (entry10_->delyavg - entry00_->delyavg);
792 
793  delysig_ = entry00_->delysig + adcota_ * (entry01_->delysig - entry00_->delysig) +
794  adcotb_ * (entry10_->delysig - entry00_->delysig);
795 
796  mpvvav_ = entry00_->mpvvav + adcota_ * (entry01_->mpvvav - entry00_->mpvvav) +
797  adcotb_ * (entry10_->mpvvav - entry00_->mpvvav);
798 
799  sigmavav_ = entry00_->sigmavav + adcota_ * (entry01_->sigmavav - entry00_->sigmavav) +
800  adcotb_ * (entry10_->sigmavav - entry00_->sigmavav);
801 
802  kappavav_ = entry00_->kappavav + adcota_ * (entry01_->kappavav - entry00_->kappavav) +
803  adcotb_ * (entry10_->kappavav - entry00_->kappavav);
804 
805  for (int i = 0; i < 4; ++i) {
806  scalex_[i] = entry00_->scalex[i] + adcota_ * (entry01_->scalex[i] - entry00_->scalex[i]) +
807  adcotb_ * (entry10_->scalex[i] - entry00_->scalex[i]);
808 
809  scaley_[i] = entry00_->scaley[i] + adcota_ * (entry01_->scaley[i] - entry00_->scaley[i]) +
810  adcotb_ * (entry10_->scaley[i] - entry00_->scaley[i]);
811 
812  offsetx_[i] = entry00_->offsetx[i] + adcota_ * (entry01_->offsetx[i] - entry00_->offsetx[i]) +
813  adcotb_ * (entry10_->offsetx[i] - entry00_->offsetx[i]);
814  if (flip_x_)
815  offsetx_[i] = -offsetx_[i];
816 
817  offsety_[i] = entry00_->offsety[i] + adcota_ * (entry01_->offsety[i] - entry00_->offsety[i]) +
818  adcotb_ * (entry10_->offsety[i] - entry00_->offsety[i]);
819  if (flip_y_)
820  offsety_[i] = -offsety_[i];
821  }
822 
823  for (int i = 0; i < 2; ++i) {
824  for (int j = 0; j < 5; ++j) {
825  // Charge loss switches sides when cot(beta) changes sign
826  if (flip_y_) {
827  xypary0x0_[1 - i][j] = (float)entry00_->xypar[i][j];
828  xypary1x0_[1 - i][j] = (float)entry10_->xypar[i][j];
829  xypary0x1_[1 - i][j] = (float)entry01_->xypar[i][j];
830  lanpar_[1 - i][j] = entry00_->lanpar[i][j] + adcota_ * (entry01_->lanpar[i][j] - entry00_->lanpar[i][j]) +
831  adcotb_ * (entry10_->lanpar[i][j] - entry00_->lanpar[i][j]);
832  } else {
833  xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
834  xypary1x0_[i][j] = (float)entry10_->xypar[i][j];
835  xypary0x1_[i][j] = (float)entry01_->xypar[i][j];
836  lanpar_[i][j] = entry00_->lanpar[i][j] + adcota_ * (entry01_->lanpar[i][j] - entry00_->lanpar[i][j]) +
837  adcotb_ * (entry10_->lanpar[i][j] - entry00_->lanpar[i][j]);
838  }
839  }
840  }
841 
842  return success_;
843 } // interpolate
844 
845 // *************************************************************************************************************************************
851 // *************************************************************************************************************************************
852 
853 #ifdef SI_PIXEL_TEMPLATE_STANDALONE
854 void SiPixelTemplate2D::sideload(SiPixelTemplateEntry2D* entry,
855  int iDtype,
856  float locBx,
857  float locBz,
858  float lorwdy,
859  float lorwdx,
860  float q50,
861  float fbin[3],
862  float xsize,
863  float ysize,
864  float zsize) {
865  // Set class variables to the input parameters
866 
867  entry00_ = entry;
868  entry01_ = entry;
869  entry10_ = entry;
870  Dtype_ = iDtype;
871  lorywidth_ = lorwdy;
872  lorxwidth_ = lorwdx;
873  xsize_ = xsize;
874  ysize_ = ysize;
875  zsize_ = zsize;
876  s50_ = q50;
877  qscale_ = 1.f;
878  for (int i = 0; i < 3; ++i) {
879  fbin_[i] = fbin[i];
880  }
881 
882  // Set other class variables
883 
884  adcota_ = 0.f;
885  adcotb_ = 0.f;
886 
887  // Interpolate things in cot(alpha)-cot(beta)
888 
889  qavg_ = entry00_->qavg;
890 
891  pixmax_ = entry00_->pixmax;
892 
893  sxymax_ = entry00_->sxymax;
894 
895  clsleny_ = entry00_->clsleny;
896 
897  clslenx_ = entry00_->clslenx;
898 
899  scaleyavg_ = 1.f;
900 
901  scalexavg_ = 1.f;
902 
903  delyavg_ = 0.f;
904 
905  delysig_ = 0.f;
906 
907  for (int i = 0; i < 4; ++i) {
908  scalex_[i] = 1.f;
909  scaley_[i] = 1.f;
910  offsetx_[i] = 0.f;
911  offsety_[i] = 0.f;
912  }
913 
914  // This works only for IP-related tracks
915 
916  flip_x_ = false;
917  flip_y_ = false;
918  float cotbeta = entry00_->cotbeta;
919  switch (Dtype_) {
920  case 0:
921  if (cotbeta < 0.f) {
922  flip_y_ = true;
923  }
924  break;
925  case 1:
926  if (locBz > 0.f) {
927  flip_y_ = true;
928  }
929  break;
930  case 2:
931  case 3:
932  case 4:
933  case 5:
934  if (locBx * locBz < 0.f) {
935  flip_x_ = true;
936  }
937  if (locBx < 0.f) {
938  flip_y_ = true;
939  }
940  break;
941  default:
942  std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << iDtype << std::endl;
943  }
944 
945  // Calculate signed quantities
946 
947  lorydrift_ = lorywidth_ / 2.;
948  if (flip_y_)
949  lorydrift_ = -lorydrift_;
950  lorxdrift_ = lorxwidth_ / 2.;
951  if (flip_x_)
952  lorxdrift_ = -lorxdrift_;
953 
954  for (int i = 0; i < 2; ++i) {
955  for (int j = 0; j < 5; ++j) {
956  // Charge loss switches sides when cot(beta) changes sign
957  if (flip_y_) {
958  xypary0x0_[1 - i][j] = (float)entry00_->xypar[i][j];
959  xypary1x0_[1 - i][j] = (float)entry00_->xypar[i][j];
960  xypary0x1_[1 - i][j] = (float)entry00_->xypar[i][j];
961  lanpar_[1 - i][j] = entry00_->lanpar[i][j];
962  } else {
963  xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
964  xypary1x0_[i][j] = (float)entry00_->xypar[i][j];
965  xypary0x1_[i][j] = (float)entry00_->xypar[i][j];
966  lanpar_[i][j] = entry00_->lanpar[i][j];
967  }
968  }
969  }
970  return;
971 }
972 #endif
973 
974 // *************************************************************************************************************************************
980 // *************************************************************************************************************************************
981 
983  float yhit,
984  bool ydouble[BYM2],
985  bool xdouble[BXM2],
986  float template2d[BXM2][BYM2],
987  bool derivatives,
988  float dpdx2d[2][BXM2][BYM2],
989  float& QTemplate) {
990  // Interpolate for a new set of track angles
991 
992  // Local variables
993  int pixx, pixy, k0, k1, l0, l1, deltax, deltay, iflipy, jflipx, imin, imax, jmin, jmax;
994  int m, n;
995  float dx, dy, ddx, ddy, adx, ady;
996  // const float deltaxy[2] = {8.33f, 12.5f};
997  const float deltaxy[2] = {16.67f, 25.0f};
998 
999  // Check to see if interpolation is valid
1000 
1001  // next, determine the indices of the closest point in k (y-displacement), l (x-displacement)
1002  // pixy and pixx are the indices of the struck pixel in the (Ty,Tx) system
1003  // k0,k1 are the k-indices of the closest and next closest point
1004  // l0,l1 are the l-indices of the closest and next closest point
1005 
1006  pixy = (int)floorf(yhit / ysize_);
1007  dy = yhit - (pixy + 0.5f) * ysize_;
1008  if (flip_y_) {
1009  dy = -dy;
1010  }
1011  k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1012  if (k0 < 0)
1013  k0 = 0;
1014  if (k0 > 6)
1015  k0 = 6;
1016  ddy = 6.f * dy / ysize_ - (k0 - 3);
1017  ady = fabs(ddy);
1018  if (ddy > 0.f) {
1019  k1 = k0 + 1;
1020  if (k1 > 6)
1021  k1 = k0 - 1;
1022  } else {
1023  k1 = k0 - 1;
1024  if (k1 < 0)
1025  k1 = k0 + 1;
1026  }
1027  pixx = (int)floorf(xhit / xsize_);
1028  dx = xhit - (pixx + 0.5f) * xsize_;
1029  if (flip_x_) {
1030  dx = -dx;
1031  }
1032  l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1033  if (l0 < 0)
1034  l0 = 0;
1035  if (l0 > 6)
1036  l0 = 6;
1037  ddx = 6.f * dx / xsize_ - (l0 - 3);
1038  adx = fabs(ddx);
1039  if (ddx > 0.f) {
1040  l1 = l0 + 1;
1041  if (l1 > 6)
1042  l1 = l0 - 1;
1043  } else {
1044  l1 = l0 - 1;
1045  if (l1 < 0)
1046  l1 = l0 + 1;
1047  }
1048 
1049  // OK, lets do the template interpolation.
1050 
1051  // First find the limits of the indices for non-zero pixels
1052 
1053  imin = std::min(entry00_->iymin, entry10_->iymin);
1054  imin_ = std::min(imin, entry01_->iymin);
1055 
1056  jmin = std::min(entry00_->jxmin, entry10_->jxmin);
1057  jmin_ = std::min(jmin, entry01_->jxmin);
1058 
1059  imax = std::max(entry00_->iymax, entry10_->iymax);
1060  imax_ = std::max(imax, entry01_->iymax);
1061 
1062  jmax = std::max(entry00_->jxmax, entry10_->jxmax);
1063  jmax_ = std::max(jmax, entry01_->jxmax);
1064 
1065  // Calculate the x and y offsets to make the new template
1066 
1067  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1068 
1069  ++pixy;
1070  ++pixx;
1071 
1072  // In the template store, the struck pixel is always (THy,THx)
1073 
1074  deltax = pixx - T2HX;
1075  deltay = pixy - T2HY;
1076 
1077  // First zero the local 2-d template
1078 
1079  for (int j = 0; j < BXM2; ++j) {
1080  for (int i = 0; i < BYM2; ++i) {
1081  xytemp_[j][i] = 0.f;
1082  }
1083  }
1084 
1085  // Loop over the non-zero part of the template index space and interpolate
1086 
1087  for (int j = jmin_; j <= jmax_; ++j) {
1088  // Flip indices as needed
1089  if (flip_x_) {
1090  jflipx = T2XSIZE - 1 - j;
1091  m = deltax + jflipx;
1092  } else {
1093  m = deltax + j;
1094  }
1095  for (int i = imin_; i <= imax_; ++i) {
1096  if (flip_y_) {
1097  iflipy = T2YSIZE - 1 - i;
1098  n = deltay + iflipy;
1099  } else {
1100  n = deltay + i;
1101  }
1102  if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1103  xytemp_[m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1104  adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1105  ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1106  adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1107  adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1108  }
1109  }
1110  }
1111 
1112  //combine rows and columns to simulate double pixels
1113 
1114  for (int n = 1; n < BYM3; ++n) {
1115  if (ydouble[n]) {
1116  // Combine the y-columns
1117  for (int m = 1; m < BXM3; ++m) {
1118  xytemp_[m][n] += xytemp_[m][n + 1];
1119  }
1120  // Now shift the remaining pixels over by one column
1121  for (int i = n + 1; i < BYM3; ++i) {
1122  for (int m = 1; m < BXM3; ++m) {
1123  xytemp_[m][i] = xytemp_[m][i + 1];
1124  }
1125  }
1126  }
1127  }
1128 
1129  //combine rows and columns to simulate double pixels
1130 
1131  for (int m = 1; m < BXM3; ++m) {
1132  if (xdouble[m]) {
1133  // Combine the x-rows
1134  for (int n = 1; n < BYM3; ++n) {
1135  xytemp_[m][n] += xytemp_[m + 1][n];
1136  }
1137  // Now shift the remaining pixels over by one row
1138  for (int j = m + 1; j < BXM3; ++j) {
1139  for (n = 1; n < BYM3; ++n) {
1140  xytemp_[j][n] = xytemp_[j + 1][n];
1141  }
1142  }
1143  }
1144  }
1145 
1146  // Finally, loop through and increment the external template
1147 
1148  float qtemptot = 0.f;
1149 
1150  for (int n = 1; n < BYM3; ++n) {
1151  for (int m = 1; m < BXM3; ++m) {
1152  if (xytemp_[m][n] != 0.f) {
1153  template2d[m][n] += xytemp_[m][n];
1154  qtemptot += xytemp_[m][n];
1155  }
1156  }
1157  }
1158 
1159  QTemplate = qtemptot;
1160 
1161  if (derivatives) {
1162  float dxytempdx[2][BXM2][BYM2], dxytempdy[2][BXM2][BYM2];
1163 
1164  for (int k = 0; k < 2; ++k) {
1165  for (int i = 0; i < BXM2; ++i) {
1166  for (int j = 0; j < BYM2; ++j) {
1167  dxytempdx[k][i][j] = 0.f;
1168  dxytempdy[k][i][j] = 0.f;
1169  dpdx2d[k][i][j] = 0.f;
1170  }
1171  }
1172  }
1173 
1174  // First do shifted +x template
1175 
1176  pixx = (int)floorf((xhit + deltaxy[0]) / xsize_);
1177  dx = (xhit + deltaxy[0]) - (pixx + 0.5f) * xsize_;
1178  if (flip_x_) {
1179  dx = -dx;
1180  }
1181  l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1182  if (l0 < 0)
1183  l0 = 0;
1184  if (l0 > 6)
1185  l0 = 6;
1186  ddx = 6.f * dx / xsize_ - (l0 - 3);
1187  adx = fabs(ddx);
1188  if (ddx > 0.f) {
1189  l1 = l0 + 1;
1190  if (l1 > 6)
1191  l1 = l0 - 1;
1192  } else {
1193  l1 = l0 - 1;
1194  if (l1 < 0)
1195  l1 = l0 + 1;
1196  }
1197 
1198  // OK, lets do the template interpolation.
1199 
1200  // Calculate the x and y offsets to make the new template
1201 
1202  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1203 
1204  ++pixx;
1205 
1206  // In the template store, the struck pixel is always (THy,THx)
1207 
1208  deltax = pixx - T2HX;
1209 
1210  // Loop over the non-zero part of the template index space and interpolate
1211 
1212  for (int j = jmin_; j <= jmax_; ++j) {
1213  // Flip indices as needed
1214  if (flip_x_) {
1215  jflipx = T2XSIZE - 1 - j;
1216  m = deltax + jflipx;
1217  } else {
1218  m = deltax + j;
1219  }
1220  for (int i = imin_; i <= imax_; ++i) {
1221  if (flip_y_) {
1222  iflipy = T2YSIZE - 1 - i;
1223  n = deltay + iflipy;
1224  } else {
1225  n = deltay + i;
1226  }
1227  if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1228  dxytempdx[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1229  adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1230  ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1231  adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1232  adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1233  }
1234  }
1235  }
1236 
1237  //combine rows and columns to simulate double pixels
1238 
1239  for (int n = 1; n < BYM3; ++n) {
1240  if (ydouble[n]) {
1241  // Combine the y-columns
1242  for (int m = 1; m < BXM3; ++m) {
1243  dxytempdx[1][m][n] += dxytempdx[1][m][n + 1];
1244  }
1245  // Now shift the remaining pixels over by one column
1246  for (int i = n + 1; i < BYM3; ++i) {
1247  for (int m = 1; m < BXM3; ++m) {
1248  dxytempdx[1][m][i] = dxytempdx[1][m][i + 1];
1249  }
1250  }
1251  }
1252  }
1253 
1254  //combine rows and columns to simulate double pixels
1255 
1256  for (int m = 1; m < BXM3; ++m) {
1257  if (xdouble[m]) {
1258  // Combine the x-rows
1259  for (int n = 1; n < BYM3; ++n) {
1260  dxytempdx[1][m][n] += dxytempdx[1][m + 1][n];
1261  }
1262  // Now shift the remaining pixels over by one row
1263  for (int j = m + 1; j < BXM3; ++j) {
1264  for (int n = 1; n < BYM3; ++n) {
1265  dxytempdx[1][j][n] = dxytempdx[1][j + 1][n];
1266  }
1267  }
1268  }
1269  }
1270 
1271  // Next do shifted -x template
1272 
1273  pixx = (int)floorf((xhit - deltaxy[0]) / xsize_);
1274  dx = (xhit - deltaxy[0]) - (pixx + 0.5f) * xsize_;
1275  if (flip_x_) {
1276  dx = -dx;
1277  }
1278  l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1279  if (l0 < 0)
1280  l0 = 0;
1281  if (l0 > 6)
1282  l0 = 6;
1283  ddx = 6.f * dx / xsize_ - (l0 - 3);
1284  adx = fabs(ddx);
1285  if (ddx > 0.f) {
1286  l1 = l0 + 1;
1287  if (l1 > 6)
1288  l1 = l0 - 1;
1289  } else {
1290  l1 = l0 - 1;
1291  if (l1 < 0)
1292  l1 = l0 + 1;
1293  }
1294 
1295  // OK, lets do the template interpolation.
1296 
1297  // Calculate the x and y offsets to make the new template
1298 
1299  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1300 
1301  ++pixx;
1302 
1303  // In the template store, the struck pixel is always (THy,THx)
1304 
1305  deltax = pixx - T2HX;
1306 
1307  // Loop over the non-zero part of the template index space and interpolate
1308 
1309  for (int j = jmin_; j <= jmax_; ++j) {
1310  // Flip indices as needed
1311  if (flip_x_) {
1312  jflipx = T2XSIZE - 1 - j;
1313  m = deltax + jflipx;
1314  } else {
1315  m = deltax + j;
1316  }
1317  for (int i = imin_; i <= imax_; ++i) {
1318  if (flip_y_) {
1319  iflipy = T2YSIZE - 1 - i;
1320  n = deltay + iflipy;
1321  } else {
1322  n = deltay + i;
1323  }
1324  if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1325  dxytempdx[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1326  adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1327  ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1328  adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1329  adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1330  }
1331  }
1332  }
1333 
1334  //combine rows and columns to simulate double pixels
1335 
1336  for (int n = 1; n < BYM3; ++n) {
1337  if (ydouble[n]) {
1338  // Combine the y-columns
1339  for (int m = 1; m < BXM3; ++m) {
1340  dxytempdx[0][m][n] += dxytempdx[0][m][n + 1];
1341  }
1342  // Now shift the remaining pixels over by one column
1343  for (int i = n + 1; i < BYM3; ++i) {
1344  for (int m = 1; m < BXM3; ++m) {
1345  dxytempdx[0][m][i] = dxytempdx[0][m][i + 1];
1346  }
1347  }
1348  }
1349  }
1350 
1351  //combine rows and columns to simulate double pixels
1352 
1353  for (int m = 1; m < BXM3; ++m) {
1354  if (xdouble[m]) {
1355  // Combine the x-rows
1356  for (int n = 1; n < BYM3; ++n) {
1357  dxytempdx[0][m][n] += dxytempdx[0][m + 1][n];
1358  }
1359  // Now shift the remaining pixels over by one row
1360  for (int j = m + 1; j < BXM3; ++j) {
1361  for (int n = 1; n < BYM3; ++n) {
1362  dxytempdx[0][j][n] = dxytempdx[0][j + 1][n];
1363  }
1364  }
1365  }
1366  }
1367 
1368  // Finally, normalize the derivatives and copy the results to the output array
1369 
1370  for (int n = 1; n < BYM3; ++n) {
1371  for (int m = 1; m < BXM3; ++m) {
1372  dpdx2d[0][m][n] = (dxytempdx[1][m][n] - dxytempdx[0][m][n]) / (2. * deltaxy[0]);
1373  }
1374  }
1375 
1376  // Next, do shifted y template
1377 
1378  pixy = (int)floorf((yhit + deltaxy[1]) / ysize_);
1379  dy = (yhit + deltaxy[1]) - (pixy + 0.5f) * ysize_;
1380  if (flip_y_) {
1381  dy = -dy;
1382  }
1383  k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1384  if (k0 < 0)
1385  k0 = 0;
1386  if (k0 > 6)
1387  k0 = 6;
1388  ddy = 6.f * dy / ysize_ - (k0 - 3);
1389  ady = fabs(ddy);
1390  if (ddy > 0.f) {
1391  k1 = k0 + 1;
1392  if (k1 > 6)
1393  k1 = k0 - 1;
1394  } else {
1395  k1 = k0 - 1;
1396  if (k1 < 0)
1397  k1 = k0 + 1;
1398  }
1399  pixx = (int)floorf(xhit / xsize_);
1400  dx = xhit - (pixx + 0.5f) * xsize_;
1401  if (flip_x_) {
1402  dx = -dx;
1403  }
1404  l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1405  if (l0 < 0)
1406  l0 = 0;
1407  if (l0 > 6)
1408  l0 = 6;
1409  ddx = 6.f * dx / xsize_ - (l0 - 3);
1410  adx = fabs(ddx);
1411  if (ddx > 0.f) {
1412  l1 = l0 + 1;
1413  if (l1 > 6)
1414  l1 = l0 - 1;
1415  } else {
1416  l1 = l0 - 1;
1417  if (l1 < 0)
1418  l1 = l0 + 1;
1419  }
1420 
1421  // OK, lets do the template interpolation.
1422 
1423  // Calculate the x and y offsets to make the new template
1424 
1425  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1426 
1427  ++pixy;
1428  ++pixx;
1429 
1430  // In the template store, the struck pixel is always (THy,THx)
1431 
1432  deltax = pixx - T2HX;
1433  deltay = pixy - T2HY;
1434 
1435  // Loop over the non-zero part of the template index space and interpolate
1436 
1437  for (int j = jmin_; j <= jmax_; ++j) {
1438  // Flip indices as needed
1439  if (flip_x_) {
1440  jflipx = T2XSIZE - 1 - j;
1441  m = deltax + jflipx;
1442  } else {
1443  m = deltax + j;
1444  }
1445  for (int i = imin_; i <= imax_; ++i) {
1446  if (flip_y_) {
1447  iflipy = T2YSIZE - 1 - i;
1448  n = deltay + iflipy;
1449  } else {
1450  n = deltay + i;
1451  }
1452  if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1453  dxytempdy[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1454  adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1455  ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1456  adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1457  adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1458  }
1459  }
1460  }
1461 
1462  //combine rows and columns to simulate double pixels
1463 
1464  for (int n = 1; n < BYM3; ++n) {
1465  if (ydouble[n]) {
1466  // Combine the y-columns
1467  for (int m = 1; m < BXM3; ++m) {
1468  dxytempdy[1][m][n] += dxytempdy[1][m][n + 1];
1469  }
1470  // Now shift the remaining pixels over by one column
1471  for (int i = n + 1; i < BYM3; ++i) {
1472  for (int m = 1; m < BXM3; ++m) {
1473  dxytempdy[1][m][i] = dxytempdy[1][m][i + 1];
1474  }
1475  }
1476  }
1477  }
1478 
1479  //combine rows and columns to simulate double pixels
1480 
1481  for (int m = 1; m < BXM3; ++m) {
1482  if (xdouble[m]) {
1483  // Combine the x-rows
1484  for (int n = 1; n < BYM3; ++n) {
1485  dxytempdy[1][m][n] += dxytempdy[1][m + 1][n];
1486  }
1487  // Now shift the remaining pixels over by one row
1488  for (int j = m + 1; j < BXM3; ++j) {
1489  for (int n = 1; n < BYM3; ++n) {
1490  dxytempdy[1][j][n] = dxytempdy[1][j + 1][n];
1491  }
1492  }
1493  }
1494  }
1495 
1496  // Next, do shifted -y template
1497 
1498  pixy = (int)floorf((yhit - deltaxy[1]) / ysize_);
1499  dy = (yhit - deltaxy[1]) - (pixy + 0.5f) * ysize_;
1500  if (flip_y_) {
1501  dy = -dy;
1502  }
1503  k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1504  if (k0 < 0)
1505  k0 = 0;
1506  if (k0 > 6)
1507  k0 = 6;
1508  ddy = 6.f * dy / ysize_ - (k0 - 3);
1509  ady = fabs(ddy);
1510  if (ddy > 0.f) {
1511  k1 = k0 + 1;
1512  if (k1 > 6)
1513  k1 = k0 - 1;
1514  } else {
1515  k1 = k0 - 1;
1516  if (k1 < 0)
1517  k1 = k0 + 1;
1518  }
1519 
1520  // OK, lets do the template interpolation.
1521 
1522  // Calculate the x and y offsets to make the new template
1523 
1524  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1525 
1526  ++pixy;
1527 
1528  // In the template store, the struck pixel is always (THy,THx)
1529 
1530  deltay = pixy - T2HY;
1531 
1532  // Loop over the non-zero part of the template index space and interpolate
1533 
1534  for (int j = jmin_; j <= jmax_; ++j) {
1535  // Flip indices as needed
1536  if (flip_x_) {
1537  jflipx = T2XSIZE - 1 - j;
1538  m = deltax + jflipx;
1539  } else {
1540  m = deltax + j;
1541  }
1542  for (int i = imin_; i <= imax_; ++i) {
1543  if (flip_y_) {
1544  iflipy = T2YSIZE - 1 - i;
1545  n = deltay + iflipy;
1546  } else {
1547  n = deltay + i;
1548  }
1549  if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1550  dxytempdy[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1551  adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1552  ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1553  adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1554  adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1555  }
1556  }
1557  }
1558 
1559  //combine rows and columns to simulate double pixels
1560 
1561  for (int n = 1; n < BYM3; ++n) {
1562  if (ydouble[n]) {
1563  // Combine the y-columns
1564  for (int m = 1; m < BXM3; ++m) {
1565  dxytempdy[0][m][n] += dxytempdy[0][m][n + 1];
1566  }
1567  // Now shift the remaining pixels over by one column
1568  for (int i = n + 1; i < BYM3; ++i) {
1569  for (int m = 1; m < BXM3; ++m) {
1570  dxytempdy[0][m][i] = dxytempdy[0][m][i + 1];
1571  }
1572  }
1573  }
1574  }
1575 
1576  //combine rows and columns to simulate double pixels
1577 
1578  for (int m = 1; m < BXM3; ++m) {
1579  if (xdouble[m]) {
1580  // Combine the x-rows
1581  for (int n = 1; n < BYM3; ++n) {
1582  dxytempdy[0][m][n] += dxytempdy[0][m + 1][n];
1583  }
1584  // Now shift the remaining pixels over by one row
1585  for (int j = m + 1; j < BXM3; ++j) {
1586  for (int n = 1; n < BYM3; ++n) {
1587  dxytempdy[0][j][n] = dxytempdy[0][j + 1][n];
1588  }
1589  }
1590  }
1591  }
1592 
1593  // Finally, normalize the derivatives and copy the results to the output array
1594 
1595  for (int n = 1; n < BYM3; ++n) {
1596  for (int m = 1; m < BXM3; ++m) {
1597  dpdx2d[1][m][n] = (dxytempdy[1][m][n] - dxytempdy[0][m][n]) / (2. * deltaxy[1]);
1598  }
1599  }
1600  }
1601 
1602  return success_;
1603 } // xytemp
1604 
1605 // *************************************************************************************************************************************
1612 // *************************************************************************************************************************************
1613 
1615  float xhit, float yhit, bool ydouble[BYM2], bool xdouble[BXM2], float template2d[BXM2][BYM2]) {
1616  // Interpolate for a new set of track angles
1617 
1618  bool derivatives = false;
1619  float dpdx2d[2][BXM2][BYM2];
1620  float QTemplate;
1621 
1622  return SiPixelTemplate2D::xytemp(xhit, yhit, ydouble, xdouble, template2d, derivatives, dpdx2d, QTemplate);
1623 
1624 } // xytemp
1625 
1626 // *************************************************************************************************************************************
1636 // *************************************************************************************************************************************
1637 
1639  float cotalpha,
1640  float cotbeta,
1641  float xhit,
1642  float yhit,
1643  std::vector<bool>& ydouble,
1644  std::vector<bool>& xdouble,
1645  float template2d[BXM2][BYM2]) {
1646  // Local variables
1647 
1648  bool derivatives = false;
1649  float dpdx2d[2][BXM2][BYM2];
1650  float QTemplate;
1651  float locBx = 1.f;
1652  if (cotbeta < 0.f) {
1653  locBx = -1.f;
1654  }
1655  float locBz = locBx;
1656  if (cotalpha < 0.f) {
1657  locBz = -locBx;
1658  }
1659 
1660  bool yd[BYM2], xd[BXM2];
1661 
1662  yd[0] = false;
1663  yd[BYM2 - 1] = false;
1664  for (int i = 0; i < TYSIZE; ++i) {
1665  yd[i + 1] = ydouble[i];
1666  }
1667  xd[0] = false;
1668  xd[BXM2 - 1] = false;
1669  for (int j = 0; j < TXSIZE; ++j) {
1670  xd[j + 1] = xdouble[j];
1671  }
1672 
1673  // Interpolate for a new set of track angles
1674 
1675  if (SiPixelTemplate2D::interpolate(id, cotalpha, cotbeta, locBz, locBx)) {
1676  return SiPixelTemplate2D::xytemp(xhit, yhit, yd, xd, template2d, derivatives, dpdx2d, QTemplate);
1677  } else {
1678  return false;
1679  }
1680 
1681 } // xytemp
1682 
1683 // ************************************************************************************************************
1689 // ************************************************************************************************************
1690 void SiPixelTemplate2D::xysigma2(float qpixel, int index, float& xysig2)
1691 
1692 {
1693  // Interpolate using quantities already stored in the private variables
1694 
1695  // Local variables
1696  float sigi, sigi2, sigi3, sigi4, qscale, err2, err00;
1697 
1698  // Make sure that input is OK
1699 
1700 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1701  if (index < 1 || index >= BYM2) {
1702  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::ysigma2 called with index = " << index << std::endl;
1703  }
1704 #else
1705  assert(index > 0 && index < BYM2);
1706 #endif
1707 
1708  // Define the maximum signal to use in the parameterization
1709 
1710  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1711 
1712  if (qpixel < sxymax_) {
1713  sigi = qpixel;
1714  qscale = 1.f;
1715  } else {
1716  sigi = sxymax_;
1717  qscale = qpixel / sxymax_;
1718  }
1719  sigi2 = sigi * sigi;
1720  sigi3 = sigi2 * sigi;
1721  sigi4 = sigi3 * sigi;
1722  if (index <= T2HYP1) {
1723  err00 = xypary0x0_[0][0] + xypary0x0_[0][1] * sigi + xypary0x0_[0][2] * sigi2 + xypary0x0_[0][3] * sigi3 +
1724  xypary0x0_[0][4] * sigi4;
1725  err2 = err00 +
1726  adcota_ * (xypary0x1_[0][0] + xypary0x1_[0][1] * sigi + xypary0x1_[0][2] * sigi2 + xypary0x1_[0][3] * sigi3 +
1727  xypary0x1_[0][4] * sigi4 - err00) +
1728  adcotb_ * (xypary1x0_[0][0] + xypary1x0_[0][1] * sigi + xypary1x0_[0][2] * sigi2 + xypary1x0_[0][3] * sigi3 +
1729  xypary1x0_[0][4] * sigi4 - err00);
1730  } else {
1731  err00 = xypary0x0_[1][0] + xypary0x0_[1][1] * sigi + xypary0x0_[1][2] * sigi2 + xypary0x0_[1][3] * sigi3 +
1732  xypary0x0_[1][4] * sigi4;
1733  err2 = err00 +
1734  adcota_ * (xypary0x1_[1][0] + xypary0x1_[1][1] * sigi + xypary0x1_[1][2] * sigi2 + xypary0x1_[1][3] * sigi3 +
1735  xypary0x1_[1][4] * sigi4 - err00) +
1736  adcotb_ * (xypary1x0_[1][0] + xypary1x0_[1][1] * sigi + xypary1x0_[1][2] * sigi2 + xypary1x0_[1][3] * sigi3 +
1737  xypary1x0_[1][4] * sigi4 - err00);
1738  }
1739  xysig2 = qscale * err2;
1740  if (xysig2 <= 0.f) {
1741  xysig2 = s50_ * s50_;
1742  }
1743 
1744  return;
1745 
1746 } // End xysigma2
1747 
1748 // ************************************************************************************************************
1750 // ************************************************************************************************************
1751 void SiPixelTemplate2D::landau_par(float lanpar[2][5])
1752 
1753 {
1754  // Interpolate using quantities already stored in the private variables
1755 
1756  for (int i = 0; i < 2; ++i) {
1757  for (int j = 0; j < 5; ++j) {
1758  lanpar[i][j] = lanpar_[i][j];
1759  }
1760  }
1761  return;
1762 
1763 } // End lan_par
T2HY
#define T2HY
Definition: SiPixelTemplateDefs.h:40
BXM3
#define BXM3
Definition: SiPixelTemplateDefs.h:37
SiPixelTemplateHeader2D::fbin
float fbin[3]
The QBin definitions in Q_clus/Q_avg.
Definition: SiPixelTemplate2D.h:96
SiPixelTemplate2D::xytemp
bool xytemp(float xhit, float yhit, bool ydouble[21+2], bool xdouble[13+2], float template2d[13+2][21+2], bool dervatives, float dpdx2d[2][13+2][21+2], float &QTemplate)
Definition: SiPixelTemplate2D.cc:982
mps_fire.i
i
Definition: mps_fire.py:428
SiPixel2DTemplateDBObject
Definition: SiPixel2DTemplateDBObject.h:16
T2YSIZE
#define T2YSIZE
Definition: SiPixelTemplateDefs.h:38
MessageLogger.h
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
SiPixelTemplateStore2D::entry
std::vector< std::vector< SiPixelTemplateEntry2D > > entry
Definition: SiPixelTemplate2D.h:106
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
SiPixelTemplateHeader2D::Bfield
float Bfield
Bfield in Tesla.
Definition: SiPixelTemplate2D.h:95
SiPixelTemplate2D::xysigma2
void xysigma2(float qpixel, int index, float &xysig2)
Definition: SiPixelTemplate2D.cc:1690
reco::ParticleMasses::k0
const double k0
Definition: ParticleMasses.h:7
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
SiPixelTemplateHeader2D::lorybias
float lorybias
estimate of y-lorentz bias
Definition: SiPixelTemplate2D.h:86
T2HX
#define T2HX
Definition: SiPixelTemplateDefs.h:42
min
T min(T a, T b)
Definition: MathUtil.h:58
edm
HLT enums.
Definition: AlignableModifier.h:19
mps_splice.entry
entry
Definition: mps_splice.py:68
SiPixel2DTemplateDBObject::char2float
Definition: SiPixel2DTemplateDBObject.h:70
gather_cfg.cout
cout
Definition: gather_cfg.py:144
SiPixelTemplateEntry2D
Definition: SiPixelTemplate2D.h:42
SiPixelTemplate2D::getid
bool getid(int id)
Definition: SiPixelTemplate2D.cc:535
cms::cuda::assert
assert(be >=bs)
SiPixelTemplateHeader2D::qscale
float qscale
Charge scaling to match cmssw and pixelav.
Definition: SiPixelTemplate2D.h:83
ysize
const Int_t ysize
Definition: trackSplitPlot.h:43
SiPixelTemplateHeader2D::temperature
float temperature
detector temperature in deg K
Definition: SiPixelTemplate2D.h:89
SiPixelTemplateHeader2D::ID
int ID
< template header structure
Definition: SiPixelTemplate2D.h:78
xsize
const Int_t xsize
Definition: trackSplitPlot.h:42
SiPixelTemplateHeader2D::s50
float s50
1/2 of the multihit dcol threshold in electrons
Definition: SiPixelTemplate2D.h:91
groupFilesInBlocks.temp
list temp
Definition: groupFilesInBlocks.py:142
SiPixelTemplateHeader2D::templ_version
int templ_version
Version number of the template to ensure code compatibility.
Definition: SiPixelTemplate2D.h:94
TXSIZE
#define TXSIZE
Definition: SiPixelTemplateDefs.h:30
FileInPath.h
SiPixelTemplateStore2D::head
SiPixelTemplateHeader2D head
< template storage structure
Definition: SiPixelTemplate2D.h:103
dqmiodatasetharvest.db
db
Definition: dqmiodatasetharvest.py:119
edm::FileInPath
Definition: FileInPath.h:64
SiPixelTemplateHeader2D::lorxbias
float lorxbias
estimate of x-lorentz bias
Definition: SiPixelTemplate2D.h:87
SiPixelTemplateHeader2D::ysize
float ysize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate2D.h:98
SiPixelTemplateHeader2D::lorxwidth
float lorxwidth
estimate of x-lorentz width for optimal resolution
Definition: SiPixelTemplate2D.h:85
BYM2
#define BYM2
Definition: SiPixelTemplateDefs.h:28
visualization-live-secondInstance_cfg.m
m
Definition: visualization-live-secondInstance_cfg.py:72
SiPixelTemplateHeader2D::ss50
float ss50
1/2 of the single hit dcol threshold in electrons
Definition: SiPixelTemplate2D.h:92
BXM2
#define BXM2
Definition: SiPixelTemplateDefs.h:36
SiPixelTemplateHeader2D::zsize
float zsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate2D.h:99
dqmdumpme.k
k
Definition: dqmdumpme.py:60
SiPixelTemplateHeader2D::title
char title[80]
template title
Definition: SiPixelTemplate2D.h:93
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
SiPixelTemplateHeader2D::lorywidth
float lorywidth
estimate of y-lorentz width for optimal resolution
Definition: SiPixelTemplate2D.h:84
runGCPTkAlMap.in_file
in_file
Definition: runGCPTkAlMap.py:74
SiPixelTemplate2D.h
BYM3
#define BYM3
Definition: SiPixelTemplateDefs.h:29
SiPixelTemplateStore2D
Definition: SiPixelTemplate2D.h:102
SiPixelTemplate2D::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
Definition: SiPixelTemplate2D.cc:625
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
T2XSIZE
#define T2XSIZE
Definition: SiPixelTemplateDefs.h:39
createfilelist.int
int
Definition: createfilelist.py:10
TYSIZE
#define TYSIZE
Definition: SiPixelTemplateDefs.h:21
FrontierConditions_GlobalTag_cff.file
file
Definition: FrontierConditions_GlobalTag_cff.py:13
SiPixelTemplateHeader2D::NTyx
int NTyx
number of Template y-slices of x entries
Definition: SiPixelTemplate2D.h:80
PVValHelper::dy
Definition: PVValidationHelpers.h:50
B2GTnPMonitor_cfi.item
item
Definition: B2GTnPMonitor_cfi.py:147
cmsLHEtoEOSManager.l
l
Definition: cmsLHEtoEOSManager.py:204
SiPixelTemplateHeader2D::NTxx
int NTxx
number of Template x-entries in each slice
Definition: SiPixelTemplate2D.h:81
SiPixelTemplateHeader2D::Dtype
int Dtype
detector type (0=BPix, 1=FPix)
Definition: SiPixelTemplate2D.h:82
SiPixelTemplate2D::landau_par
void landau_par(float lanpar[2][5])
Return the Landau probability parameters for this set of cot(alpha, cot(beta)
Definition: SiPixelTemplate2D.cc:1751
ENDL
#define ENDL
Definition: SiPixelTemplate2D.cc:45
triggerObjects_cff.id
id
Definition: triggerObjects_cff.py:29
Exception
Definition: hltDiff.cc:245
Exception.h
LOGINFO
#define LOGINFO(x)
Definition: SiPixelTemplate2D.cc:44
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
SiPixelTemplateHeader2D::NTy
int NTy
number of Template y entries
Definition: SiPixelTemplate2D.h:79
dummy
Definition: DummySelector.h:38
SiPixelTemplateHeader2D::Vbias
float Vbias
detector bias potential in Volts
Definition: SiPixelTemplate2D.h:88
c
auto & c
Definition: CAHitNtupletGeneratorKernelsImpl.h:46
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
LOGERROR
#define LOGERROR(x)
Definition: SiPixelTemplate2D.cc:43
SiPixelTemplateHeader2D::xsize
float xsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate2D.h:97
SiPixelTemplate2D::pushfile
static bool pushfile(int filenum, std::vector< SiPixelTemplateStore2D > &pixelTemp, std::string dir="CalibTracker/SiPixelESProducers/data/")
Definition: SiPixelTemplate2D.cc:61
PVValHelper::dx
Definition: PVValidationHelpers.h:49
SiPixelTemplateHeader2D::fluence
float fluence
radiation fluence in n_eq/cm^2
Definition: SiPixelTemplate2D.h:90
DeadROC_duringRun.dir
dir
Definition: DeadROC_duringRun.py:23
T2HYP1
#define T2HYP1
Definition: SiPixelTemplateDefs.h:41