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