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