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 charge scaling factor 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  success_ = getid(id);
672  }
673 
674 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
675  if (index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
676  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
677  << ", Are you using the correct global tag?" << std::endl;
678  }
679 #else
680  assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
681 #endif
682 
683  // Check angle limits and et up interpolation parameters
684 
685  float cota = cotalpha;
686  flip_x_ = false;
687  flip_y_ = false;
688  switch (Dtype_) {
689  case 0:
690  if (cotbeta < 0.f) {
691  flip_y_ = true;
692  }
693  break;
694  case 1:
695  if (locBz > 0.f) {
696  flip_y_ = true;
697  }
698  break;
699  case 2:
700  case 3:
701  case 4:
702  case 5:
703  if (locBx * locBz < 0.f) {
704  cota = fabs(cotalpha);
705  flip_x_ = true;
706  }
707  if (locBx < 0.f) {
708  flip_y_ = true;
709  }
710  break;
711  default:
712 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
713  throw cms::Exception("DataCorrupt")
714  << "SiPixelTemplate2D::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
715 #else
716  std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
717 #endif
718  }
719 
720  if (cota < cotalpha0_) {
721  success_ = false;
722  jx0_ = 0;
723  jx1_ = 1;
724  adcota_ = 0.f;
725  } else if (cota > cotalpha1_) {
726  success_ = false;
727  jx0_ = Nxx_ - 1;
728  jx1_ = jx0_ - 1;
729  adcota_ = 0.f;
730  } else {
731  jx0_ = (int)((cota - cotalpha0_) / deltacota_ + 0.5f);
732  dcota = (cota - (cotalpha0_ + jx0_ * deltacota_)) / deltacota_;
733  adcota_ = fabs(dcota);
734  if (dcota > 0.f) {
735  jx1_ = jx0_ + 1;
736  if (jx1_ > Nxx_ - 1)
737  jx1_ = jx0_ - 1;
738  } else {
739  jx1_ = jx0_ - 1;
740  if (jx1_ < 0)
741  jx1_ = jx0_ + 1;
742  }
743  }
744 
745  // Interpolate the absolute value of cot(beta)
746 
747  acotb = fabs(cotbeta);
748 
749  if (acotb < cotbeta0_) {
750  success_ = false;
751  iy0_ = 0;
752  iy1_ = 1;
753  adcotb_ = 0.f;
754  } else if (acotb > cotbeta1_) {
755  success_ = false;
756  iy0_ = Nyx_ - 1;
757  iy1_ = iy0_ - 1;
758  adcotb_ = 0.f;
759  } else {
760  iy0_ = (int)((acotb - cotbeta0_) / deltacotb_ + 0.5f);
761  dcotb = (acotb - (cotbeta0_ + iy0_ * deltacotb_)) / deltacotb_;
762  adcotb_ = fabs(dcotb);
763  if (dcotb > 0.f) {
764  iy1_ = iy0_ + 1;
765  if (iy1_ > Nyx_ - 1)
766  iy1_ = iy0_ - 1;
767  } else {
768  iy1_ = iy0_ - 1;
769  if (iy1_ < 0)
770  iy1_ = iy0_ + 1;
771  }
772  }
773 
774  // Calculate signed quantities
775 
776  lorydrift_ = lorywidth_ / 2.;
777  if (flip_y_)
778  lorydrift_ = -lorydrift_;
779  lorxdrift_ = lorxwidth_ / 2.;
780  if (flip_x_)
781  lorxdrift_ = -lorxdrift_;
782 
783  // Use pointers to the three angle pairs used in the interpolation
784 
785  entry00_ = &thePixelTemp_[index_id_].entry[iy0_][jx0_];
786  entry10_ = &thePixelTemp_[index_id_].entry[iy1_][jx0_];
787  entry01_ = &thePixelTemp_[index_id_].entry[iy0_][jx1_];
788 
789  // Interpolate things in cot(alpha)-cot(beta)
790 
791  qavg_ = entry00_->qavg + adcota_ * (entry01_->qavg - entry00_->qavg) + adcotb_ * (entry10_->qavg - entry00_->qavg);
792 
793  pixmax_ = entry00_->pixmax + adcota_ * (entry01_->pixmax - entry00_->pixmax) +
794  adcotb_ * (entry10_->pixmax - entry00_->pixmax);
795 
796  sxymax_ = entry00_->sxymax + adcota_ * (entry01_->sxymax - entry00_->sxymax) +
797  adcotb_ * (entry10_->sxymax - entry00_->sxymax);
798 
799  chi2avgone_ = entry00_->chi2avgone + adcota_ * (entry01_->chi2avgone - entry00_->chi2avgone) +
800  adcotb_ * (entry10_->chi2avgone - entry00_->chi2avgone);
801 
802  chi2minone_ = entry00_->chi2minone + adcota_ * (entry01_->chi2minone - entry00_->chi2minone) +
803  adcotb_ * (entry10_->chi2minone - entry00_->chi2minone);
804 
805  clsleny_ = entry00_->clsleny + adcota_ * (entry01_->clsleny - entry00_->clsleny) +
806  adcotb_ * (entry10_->clsleny - entry00_->clsleny);
807 
808  clslenx_ = entry00_->clslenx + adcota_ * (entry01_->clslenx - entry00_->clslenx) +
809  adcotb_ * (entry10_->clslenx - entry00_->clslenx);
810 
811  chi2ppix_ = entry00_->chi2ppix + adcota_ * (entry01_->chi2ppix - entry00_->chi2ppix) +
812  adcotb_ * (entry10_->chi2ppix - entry00_->chi2ppix);
813 
814  chi2scale_ = entry00_->chi2scale + adcota_ * (entry01_->chi2scale - entry00_->chi2scale) +
815  adcotb_ * (entry10_->chi2scale - entry00_->chi2scale);
816 
817  scaleyavg_ = entry00_->scaleyavg + adcota_ * (entry01_->scaleyavg - entry00_->scaleyavg) +
818  adcotb_ * (entry10_->scaleyavg - entry00_->scaleyavg);
819 
820  scalexavg_ = entry00_->scalexavg + adcota_ * (entry01_->scalexavg - entry00_->scalexavg) +
821  adcotb_ * (entry10_->scalexavg - entry00_->scalexavg);
822 
823  delyavg_ = entry00_->delyavg + adcota_ * (entry01_->delyavg - entry00_->delyavg) +
824  adcotb_ * (entry10_->delyavg - entry00_->delyavg);
825 
826  delysig_ = entry00_->delysig + adcota_ * (entry01_->delysig - entry00_->delysig) +
827  adcotb_ * (entry10_->delysig - entry00_->delysig);
828 
829  mpvvav_ = entry00_->mpvvav + adcota_ * (entry01_->mpvvav - entry00_->mpvvav) +
830  adcotb_ * (entry10_->mpvvav - entry00_->mpvvav);
831 
832  sigmavav_ = entry00_->sigmavav + adcota_ * (entry01_->sigmavav - entry00_->sigmavav) +
833  adcotb_ * (entry10_->sigmavav - entry00_->sigmavav);
834 
835  kappavav_ = entry00_->kappavav + adcota_ * (entry01_->kappavav - entry00_->kappavav) +
836  adcotb_ * (entry10_->kappavav - entry00_->kappavav);
837 
838  for (int i = 0; i < 4; ++i) {
839  scalex_[i] = entry00_->scalex[i] + adcota_ * (entry01_->scalex[i] - entry00_->scalex[i]) +
840  adcotb_ * (entry10_->scalex[i] - entry00_->scalex[i]);
841 
842  scaley_[i] = entry00_->scaley[i] + adcota_ * (entry01_->scaley[i] - entry00_->scaley[i]) +
843  adcotb_ * (entry10_->scaley[i] - entry00_->scaley[i]);
844 
845  offsetx_[i] = entry00_->offsetx[i] + adcota_ * (entry01_->offsetx[i] - entry00_->offsetx[i]) +
846  adcotb_ * (entry10_->offsetx[i] - entry00_->offsetx[i]);
847  if (flip_x_)
848  offsetx_[i] = -offsetx_[i];
849 
850  offsety_[i] = entry00_->offsety[i] + adcota_ * (entry01_->offsety[i] - entry00_->offsety[i]) +
851  adcotb_ * (entry10_->offsety[i] - entry00_->offsety[i]);
852  if (flip_y_)
853  offsety_[i] = -offsety_[i];
854  }
855 
856  for (int i = 0; i < 2; ++i) {
857  for (int j = 0; j < 5; ++j) {
858  // Charge loss switches sides when cot(beta) changes sign
859  if (flip_y_) {
860  xypary0x0_[1 - i][j] = (float)entry00_->xypar[i][j];
861  xypary1x0_[1 - i][j] = (float)entry10_->xypar[i][j];
862  xypary0x1_[1 - i][j] = (float)entry01_->xypar[i][j];
863  lanpar_[1 - i][j] = entry00_->lanpar[i][j] + adcota_ * (entry01_->lanpar[i][j] - entry00_->lanpar[i][j]) +
864  adcotb_ * (entry10_->lanpar[i][j] - entry00_->lanpar[i][j]);
865  } else {
866  xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
867  xypary1x0_[i][j] = (float)entry10_->xypar[i][j];
868  xypary0x1_[i][j] = (float)entry01_->xypar[i][j];
869  lanpar_[i][j] = entry00_->lanpar[i][j] + adcota_ * (entry01_->lanpar[i][j] - entry00_->lanpar[i][j]) +
870  adcotb_ * (entry10_->lanpar[i][j] - entry00_->lanpar[i][j]);
871  }
872  }
873  }
874 
875  return success_;
876 } // interpolate
877 
878 // *************************************************************************************************************************************
884 // *************************************************************************************************************************************
885 
886 void SiPixelTemplate2D::sideload(SiPixelTemplateEntry2D* entry,
887  int iDtype,
888  float locBx,
889  float locBz,
890  float lorwdy,
891  float lorwdx,
892  float q50,
893  float fbin[3],
894  float xsize,
895  float ysize,
896  float zsize) {
897  // Set class variables to the input parameters
898 
899  entry00_ = entry;
900  entry01_ = entry;
901  entry10_ = entry;
902  Dtype_ = iDtype;
903  lorywidth_ = lorwdy;
904  lorxwidth_ = lorwdx;
905  xsize_ = xsize;
906  ysize_ = ysize;
907  zsize_ = zsize;
908  s50_ = q50;
909  qscale_ = 1.f;
910  for (int i = 0; i < 3; ++i) {
911  fbin_[i] = fbin[i];
912  }
913 
914  // Set other class variables
915 
916  adcota_ = 0.f;
917  adcotb_ = 0.f;
918 
919  // Interpolate things in cot(alpha)-cot(beta)
920 
921  qavg_ = entry00_->qavg;
922 
923  pixmax_ = entry00_->pixmax;
924 
925  sxymax_ = entry00_->sxymax;
926 
927  clsleny_ = entry00_->clsleny;
928 
929  clslenx_ = entry00_->clslenx;
930 
931  scaleyavg_ = 1.f;
932 
933  scalexavg_ = 1.f;
934 
935  delyavg_ = 0.f;
936 
937  delysig_ = 0.f;
938 
939  for (int i = 0; i < 4; ++i) {
940  scalex_[i] = 1.f;
941  scaley_[i] = 1.f;
942  offsetx_[i] = 0.f;
943  offsety_[i] = 0.f;
944  }
945 
946  // This works only for IP-related tracks
947 
948  flip_x_ = false;
949  flip_y_ = false;
950  float cotbeta = entry00_->cotbeta;
951  switch (Dtype_) {
952  case 0:
953  if (cotbeta < 0.f) {
954  flip_y_ = true;
955  }
956  break;
957  case 1:
958  if (locBz > 0.f) {
959  flip_y_ = true;
960  }
961  break;
962  case 2:
963  case 3:
964  case 4:
965  case 5:
966  if (locBx * locBz < 0.f) {
967  flip_x_ = true;
968  }
969  if (locBx < 0.f) {
970  flip_y_ = true;
971  }
972  break;
973  default:
974 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
975  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::illegal subdetector ID = " << iDtype << std::endl;
976 #else
977  std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << iDtype << std::endl;
978 #endif
979  }
980 
981  // Calculate signed quantities
982 
983  lorydrift_ = lorywidth_ / 2.;
984  if (flip_y_)
985  lorydrift_ = -lorydrift_;
986  lorxdrift_ = lorxwidth_ / 2.;
987  if (flip_x_)
988  lorxdrift_ = -lorxdrift_;
989 
990  for (int i = 0; i < 2; ++i) {
991  for (int j = 0; j < 5; ++j) {
992  // Charge loss switches sides when cot(beta) changes sign
993  if (flip_y_) {
994  xypary0x0_[1 - i][j] = (float)entry00_->xypar[i][j];
995  xypary1x0_[1 - i][j] = (float)entry00_->xypar[i][j];
996  xypary0x1_[1 - i][j] = (float)entry00_->xypar[i][j];
997  lanpar_[1 - i][j] = entry00_->lanpar[i][j];
998  } else {
999  xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
1000  xypary1x0_[i][j] = (float)entry00_->xypar[i][j];
1001  xypary0x1_[i][j] = (float)entry00_->xypar[i][j];
1002  lanpar_[i][j] = entry00_->lanpar[i][j];
1003  }
1004  }
1005  }
1006  return;
1007 }
1008 
1009 // *************************************************************************************************************************************
1015 // *************************************************************************************************************************************
1016 
1017 bool SiPixelTemplate2D::xytemp(float xhit,
1018  float yhit,
1019  bool ydouble[BYM2],
1020  bool xdouble[BXM2],
1021  float template2d[BXM2][BYM2],
1022  bool derivatives,
1023  float dpdx2d[2][BXM2][BYM2],
1024  float& QTemplate) {
1025  // Interpolate for a new set of track angles
1026 
1027  // Local variables
1028  int pixx, pixy, k0, k1, l0, l1, deltax, deltay, iflipy, jflipx, imin, imax, jmin, jmax;
1029  int m, n;
1030  float dx, dy, ddx, ddy, adx, ady;
1031  // const float deltaxy[2] = {8.33f, 12.5f};
1032  const float deltaxy[2] = {16.67f, 25.0f};
1033 
1034  // Check to see if interpolation is valid
1035 
1036  // next, determine the indices of the closest point in k (y-displacement), l (x-displacement)
1037  // pixy and pixx are the indices of the struck pixel in the (Ty,Tx) system
1038  // k0,k1 are the k-indices of the closest and next closest point
1039  // l0,l1 are the l-indices of the closest and next closest point
1040 
1041  pixy = (int)floorf(yhit / ysize_);
1042  dy = yhit - (pixy + 0.5f) * ysize_;
1043  if (flip_y_) {
1044  dy = -dy;
1045  }
1046  k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1047  if (k0 < 0)
1048  k0 = 0;
1049  if (k0 > 6)
1050  k0 = 6;
1051  ddy = 6.f * dy / ysize_ - (k0 - 3);
1052  ady = fabs(ddy);
1053  if (ddy > 0.f) {
1054  k1 = k0 + 1;
1055  if (k1 > 6)
1056  k1 = k0 - 1;
1057  } else {
1058  k1 = k0 - 1;
1059  if (k1 < 0)
1060  k1 = k0 + 1;
1061  }
1062  pixx = (int)floorf(xhit / xsize_);
1063  dx = xhit - (pixx + 0.5f) * xsize_;
1064  if (flip_x_) {
1065  dx = -dx;
1066  }
1067  l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1068  if (l0 < 0)
1069  l0 = 0;
1070  if (l0 > 6)
1071  l0 = 6;
1072  ddx = 6.f * dx / xsize_ - (l0 - 3);
1073  adx = fabs(ddx);
1074  if (ddx > 0.f) {
1075  l1 = l0 + 1;
1076  if (l1 > 6)
1077  l1 = l0 - 1;
1078  } else {
1079  l1 = l0 - 1;
1080  if (l1 < 0)
1081  l1 = l0 + 1;
1082  }
1083 
1084  // OK, lets do the template interpolation.
1085 
1086  // First find the limits of the indices for non-zero pixels
1087 
1088  imin = std::min(entry00_->iymin, entry10_->iymin);
1089  imin_ = std::min(imin, entry01_->iymin);
1090 
1091  jmin = std::min(entry00_->jxmin, entry10_->jxmin);
1092  jmin_ = std::min(jmin, entry01_->jxmin);
1093 
1094  imax = std::max(entry00_->iymax, entry10_->iymax);
1095  imax_ = std::max(imax, entry01_->iymax);
1096 
1097  jmax = std::max(entry00_->jxmax, entry10_->jxmax);
1098  jmax_ = std::max(jmax, entry01_->jxmax);
1099 
1100  // Calculate the x and y offsets to make the new template
1101 
1102  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1103 
1104  ++pixy;
1105  ++pixx;
1106 
1107  // In the template store, the struck pixel is always (THy,THx)
1108 
1109  deltax = pixx - T2HX;
1110  deltay = pixy - T2HY;
1111 
1112  // First zero the local 2-d template
1113 
1114  for (int j = 0; j < BXM2; ++j) {
1115  for (int i = 0; i < BYM2; ++i) {
1116  xytemp_[j][i] = 0.f;
1117  }
1118  }
1119 
1120  // Loop over the non-zero part of the template index space and interpolate
1121 
1122  for (int j = jmin_; j <= jmax_; ++j) {
1123  // Flip indices as needed
1124  if (flip_x_) {
1125  jflipx = T2XSIZE - 1 - j;
1126  m = deltax + jflipx;
1127  } else {
1128  m = deltax + j;
1129  }
1130  for (int i = imin_; i <= imax_; ++i) {
1131  if (flip_y_) {
1132  iflipy = T2YSIZE - 1 - i;
1133  n = deltay + iflipy;
1134  } else {
1135  n = deltay + i;
1136  }
1137  if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1138  xytemp_[m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1139  adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1140  ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1141  adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1142  adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1143  }
1144  }
1145  }
1146 
1147  //combine rows and columns to simulate double pixels
1148 
1149  for (int n = 1; n < BYM3; ++n) {
1150  if (ydouble[n]) {
1151  // Combine the y-columns
1152  for (int m = 1; m < BXM3; ++m) {
1153  xytemp_[m][n] += xytemp_[m][n + 1];
1154  }
1155  // Now shift the remaining pixels over by one column
1156  for (int i = n + 1; i < BYM3; ++i) {
1157  for (int m = 1; m < BXM3; ++m) {
1158  xytemp_[m][i] = xytemp_[m][i + 1];
1159  }
1160  }
1161  }
1162  }
1163 
1164  //combine rows and columns to simulate double pixels
1165 
1166  for (int m = 1; m < BXM3; ++m) {
1167  if (xdouble[m]) {
1168  // Combine the x-rows
1169  for (int n = 1; n < BYM3; ++n) {
1170  xytemp_[m][n] += xytemp_[m + 1][n];
1171  }
1172  // Now shift the remaining pixels over by one row
1173  for (int j = m + 1; j < BXM3; ++j) {
1174  for (n = 1; n < BYM3; ++n) {
1175  xytemp_[j][n] = xytemp_[j + 1][n];
1176  }
1177  }
1178  }
1179  }
1180 
1181  // Finally, loop through and increment the external template
1182 
1183  float qtemptot = 0.f;
1184 
1185  for (int n = 1; n < BYM3; ++n) {
1186  for (int m = 1; m < BXM3; ++m) {
1187  if (xytemp_[m][n] != 0.f) {
1188  template2d[m][n] += xytemp_[m][n];
1189  qtemptot += xytemp_[m][n];
1190  }
1191  }
1192  }
1193 
1194  QTemplate = qtemptot;
1195 
1196  if (derivatives) {
1197  float dxytempdx[2][BXM2][BYM2], dxytempdy[2][BXM2][BYM2];
1198 
1199  for (int k = 0; k < 2; ++k) {
1200  for (int i = 0; i < BXM2; ++i) {
1201  for (int j = 0; j < BYM2; ++j) {
1202  dxytempdx[k][i][j] = 0.f;
1203  dxytempdy[k][i][j] = 0.f;
1204  dpdx2d[k][i][j] = 0.f;
1205  }
1206  }
1207  }
1208 
1209  // First do shifted +x template
1210 
1211  pixx = (int)floorf((xhit + deltaxy[0]) / xsize_);
1212  dx = (xhit + deltaxy[0]) - (pixx + 0.5f) * xsize_;
1213  if (flip_x_) {
1214  dx = -dx;
1215  }
1216  l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1217  if (l0 < 0)
1218  l0 = 0;
1219  if (l0 > 6)
1220  l0 = 6;
1221  ddx = 6.f * dx / xsize_ - (l0 - 3);
1222  adx = fabs(ddx);
1223  if (ddx > 0.f) {
1224  l1 = l0 + 1;
1225  if (l1 > 6)
1226  l1 = l0 - 1;
1227  } else {
1228  l1 = l0 - 1;
1229  if (l1 < 0)
1230  l1 = l0 + 1;
1231  }
1232 
1233  // OK, lets do the template interpolation.
1234 
1235  // Calculate the x and y offsets to make the new template
1236 
1237  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1238 
1239  ++pixx;
1240 
1241  // In the template store, the struck pixel is always (THy,THx)
1242 
1243  deltax = pixx - T2HX;
1244 
1245  // Loop over the non-zero part of the template index space and interpolate
1246 
1247  for (int j = jmin_; j <= jmax_; ++j) {
1248  // Flip indices as needed
1249  if (flip_x_) {
1250  jflipx = T2XSIZE - 1 - j;
1251  m = deltax + jflipx;
1252  } else {
1253  m = deltax + j;
1254  }
1255  for (int i = imin_; i <= imax_; ++i) {
1256  if (flip_y_) {
1257  iflipy = T2YSIZE - 1 - i;
1258  n = deltay + iflipy;
1259  } else {
1260  n = deltay + i;
1261  }
1262  if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1263  dxytempdx[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1264  adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1265  ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1266  adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1267  adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1268  }
1269  }
1270  }
1271 
1272  //combine rows and columns to simulate double pixels
1273 
1274  for (int n = 1; n < BYM3; ++n) {
1275  if (ydouble[n]) {
1276  // Combine the y-columns
1277  for (int m = 1; m < BXM3; ++m) {
1278  dxytempdx[1][m][n] += dxytempdx[1][m][n + 1];
1279  }
1280  // Now shift the remaining pixels over by one column
1281  for (int i = n + 1; i < BYM3; ++i) {
1282  for (int m = 1; m < BXM3; ++m) {
1283  dxytempdx[1][m][i] = dxytempdx[1][m][i + 1];
1284  }
1285  }
1286  }
1287  }
1288 
1289  //combine rows and columns to simulate double pixels
1290 
1291  for (int m = 1; m < BXM3; ++m) {
1292  if (xdouble[m]) {
1293  // Combine the x-rows
1294  for (int n = 1; n < BYM3; ++n) {
1295  dxytempdx[1][m][n] += dxytempdx[1][m + 1][n];
1296  }
1297  // Now shift the remaining pixels over by one row
1298  for (int j = m + 1; j < BXM3; ++j) {
1299  for (int n = 1; n < BYM3; ++n) {
1300  dxytempdx[1][j][n] = dxytempdx[1][j + 1][n];
1301  }
1302  }
1303  }
1304  }
1305 
1306  // Next do shifted -x template
1307 
1308  pixx = (int)floorf((xhit - deltaxy[0]) / xsize_);
1309  dx = (xhit - deltaxy[0]) - (pixx + 0.5f) * xsize_;
1310  if (flip_x_) {
1311  dx = -dx;
1312  }
1313  l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1314  if (l0 < 0)
1315  l0 = 0;
1316  if (l0 > 6)
1317  l0 = 6;
1318  ddx = 6.f * dx / xsize_ - (l0 - 3);
1319  adx = fabs(ddx);
1320  if (ddx > 0.f) {
1321  l1 = l0 + 1;
1322  if (l1 > 6)
1323  l1 = l0 - 1;
1324  } else {
1325  l1 = l0 - 1;
1326  if (l1 < 0)
1327  l1 = l0 + 1;
1328  }
1329 
1330  // OK, lets do the template interpolation.
1331 
1332  // Calculate the x and y offsets to make the new template
1333 
1334  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1335 
1336  ++pixx;
1337 
1338  // In the template store, the struck pixel is always (THy,THx)
1339 
1340  deltax = pixx - T2HX;
1341 
1342  // Loop over the non-zero part of the template index space and interpolate
1343 
1344  for (int j = jmin_; j <= jmax_; ++j) {
1345  // Flip indices as needed
1346  if (flip_x_) {
1347  jflipx = T2XSIZE - 1 - j;
1348  m = deltax + jflipx;
1349  } else {
1350  m = deltax + j;
1351  }
1352  for (int i = imin_; i <= imax_; ++i) {
1353  if (flip_y_) {
1354  iflipy = T2YSIZE - 1 - i;
1355  n = deltay + iflipy;
1356  } else {
1357  n = deltay + i;
1358  }
1359  if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1360  dxytempdx[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1361  adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1362  ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1363  adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1364  adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1365  }
1366  }
1367  }
1368 
1369  //combine rows and columns to simulate double pixels
1370 
1371  for (int n = 1; n < BYM3; ++n) {
1372  if (ydouble[n]) {
1373  // Combine the y-columns
1374  for (int m = 1; m < BXM3; ++m) {
1375  dxytempdx[0][m][n] += dxytempdx[0][m][n + 1];
1376  }
1377  // Now shift the remaining pixels over by one column
1378  for (int i = n + 1; i < BYM3; ++i) {
1379  for (int m = 1; m < BXM3; ++m) {
1380  dxytempdx[0][m][i] = dxytempdx[0][m][i + 1];
1381  }
1382  }
1383  }
1384  }
1385 
1386  //combine rows and columns to simulate double pixels
1387 
1388  for (int m = 1; m < BXM3; ++m) {
1389  if (xdouble[m]) {
1390  // Combine the x-rows
1391  for (int n = 1; n < BYM3; ++n) {
1392  dxytempdx[0][m][n] += dxytempdx[0][m + 1][n];
1393  }
1394  // Now shift the remaining pixels over by one row
1395  for (int j = m + 1; j < BXM3; ++j) {
1396  for (int n = 1; n < BYM3; ++n) {
1397  dxytempdx[0][j][n] = dxytempdx[0][j + 1][n];
1398  }
1399  }
1400  }
1401  }
1402 
1403  // Finally, normalize the derivatives and copy the results to the output array
1404 
1405  for (int n = 1; n < BYM3; ++n) {
1406  for (int m = 1; m < BXM3; ++m) {
1407  dpdx2d[0][m][n] = (dxytempdx[1][m][n] - dxytempdx[0][m][n]) / (2. * deltaxy[0]);
1408  }
1409  }
1410 
1411  // Next, do shifted y template
1412 
1413  pixy = (int)floorf((yhit + deltaxy[1]) / ysize_);
1414  dy = (yhit + deltaxy[1]) - (pixy + 0.5f) * ysize_;
1415  if (flip_y_) {
1416  dy = -dy;
1417  }
1418  k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1419  if (k0 < 0)
1420  k0 = 0;
1421  if (k0 > 6)
1422  k0 = 6;
1423  ddy = 6.f * dy / ysize_ - (k0 - 3);
1424  ady = fabs(ddy);
1425  if (ddy > 0.f) {
1426  k1 = k0 + 1;
1427  if (k1 > 6)
1428  k1 = k0 - 1;
1429  } else {
1430  k1 = k0 - 1;
1431  if (k1 < 0)
1432  k1 = k0 + 1;
1433  }
1434  pixx = (int)floorf(xhit / xsize_);
1435  dx = xhit - (pixx + 0.5f) * xsize_;
1436  if (flip_x_) {
1437  dx = -dx;
1438  }
1439  l0 = (int)(dx / xsize_ * 6.f + 3.5f);
1440  if (l0 < 0)
1441  l0 = 0;
1442  if (l0 > 6)
1443  l0 = 6;
1444  ddx = 6.f * dx / xsize_ - (l0 - 3);
1445  adx = fabs(ddx);
1446  if (ddx > 0.f) {
1447  l1 = l0 + 1;
1448  if (l1 > 6)
1449  l1 = l0 - 1;
1450  } else {
1451  l1 = l0 - 1;
1452  if (l1 < 0)
1453  l1 = l0 + 1;
1454  }
1455 
1456  // OK, lets do the template interpolation.
1457 
1458  // Calculate the x and y offsets to make the new template
1459 
1460  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1461 
1462  ++pixy;
1463  ++pixx;
1464 
1465  // In the template store, the struck pixel is always (THy,THx)
1466 
1467  deltax = pixx - T2HX;
1468  deltay = pixy - T2HY;
1469 
1470  // Loop over the non-zero part of the template index space and interpolate
1471 
1472  for (int j = jmin_; j <= jmax_; ++j) {
1473  // Flip indices as needed
1474  if (flip_x_) {
1475  jflipx = T2XSIZE - 1 - j;
1476  m = deltax + jflipx;
1477  } else {
1478  m = deltax + j;
1479  }
1480  for (int i = imin_; i <= imax_; ++i) {
1481  if (flip_y_) {
1482  iflipy = T2YSIZE - 1 - i;
1483  n = deltay + iflipy;
1484  } else {
1485  n = deltay + i;
1486  }
1487  if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1488  dxytempdy[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1489  adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1490  ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1491  adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1492  adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1493  }
1494  }
1495  }
1496 
1497  //combine rows and columns to simulate double pixels
1498 
1499  for (int n = 1; n < BYM3; ++n) {
1500  if (ydouble[n]) {
1501  // Combine the y-columns
1502  for (int m = 1; m < BXM3; ++m) {
1503  dxytempdy[1][m][n] += dxytempdy[1][m][n + 1];
1504  }
1505  // Now shift the remaining pixels over by one column
1506  for (int i = n + 1; i < BYM3; ++i) {
1507  for (int m = 1; m < BXM3; ++m) {
1508  dxytempdy[1][m][i] = dxytempdy[1][m][i + 1];
1509  }
1510  }
1511  }
1512  }
1513 
1514  //combine rows and columns to simulate double pixels
1515 
1516  for (int m = 1; m < BXM3; ++m) {
1517  if (xdouble[m]) {
1518  // Combine the x-rows
1519  for (int n = 1; n < BYM3; ++n) {
1520  dxytempdy[1][m][n] += dxytempdy[1][m + 1][n];
1521  }
1522  // Now shift the remaining pixels over by one row
1523  for (int j = m + 1; j < BXM3; ++j) {
1524  for (int n = 1; n < BYM3; ++n) {
1525  dxytempdy[1][j][n] = dxytempdy[1][j + 1][n];
1526  }
1527  }
1528  }
1529  }
1530 
1531  // Next, do shifted -y template
1532 
1533  pixy = (int)floorf((yhit - deltaxy[1]) / ysize_);
1534  dy = (yhit - deltaxy[1]) - (pixy + 0.5f) * ysize_;
1535  if (flip_y_) {
1536  dy = -dy;
1537  }
1538  k0 = (int)(dy / ysize_ * 6.f + 3.5f);
1539  if (k0 < 0)
1540  k0 = 0;
1541  if (k0 > 6)
1542  k0 = 6;
1543  ddy = 6.f * dy / ysize_ - (k0 - 3);
1544  ady = fabs(ddy);
1545  if (ddy > 0.f) {
1546  k1 = k0 + 1;
1547  if (k1 > 6)
1548  k1 = k0 - 1;
1549  } else {
1550  k1 = k0 - 1;
1551  if (k1 < 0)
1552  k1 = k0 + 1;
1553  }
1554 
1555  // OK, lets do the template interpolation.
1556 
1557  // Calculate the x and y offsets to make the new template
1558 
1559  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1560 
1561  ++pixy;
1562 
1563  // In the template store, the struck pixel is always (THy,THx)
1564 
1565  deltay = pixy - T2HY;
1566 
1567  // Loop over the non-zero part of the template index space and interpolate
1568 
1569  for (int j = jmin_; j <= jmax_; ++j) {
1570  // Flip indices as needed
1571  if (flip_x_) {
1572  jflipx = T2XSIZE - 1 - j;
1573  m = deltax + jflipx;
1574  } else {
1575  m = deltax + j;
1576  }
1577  for (int i = imin_; i <= imax_; ++i) {
1578  if (flip_y_) {
1579  iflipy = T2YSIZE - 1 - i;
1580  n = deltay + iflipy;
1581  } else {
1582  n = deltay + i;
1583  }
1584  if (m >= 0 && m <= BXM3 && n >= 0 && n <= BYM3) {
1585  dxytempdy[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j] +
1586  adx * (float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1587  ady * (float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1588  adcota_ * (float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]) +
1589  adcotb_ * (float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1590  }
1591  }
1592  }
1593 
1594  //combine rows and columns to simulate double pixels
1595 
1596  for (int n = 1; n < BYM3; ++n) {
1597  if (ydouble[n]) {
1598  // Combine the y-columns
1599  for (int m = 1; m < BXM3; ++m) {
1600  dxytempdy[0][m][n] += dxytempdy[0][m][n + 1];
1601  }
1602  // Now shift the remaining pixels over by one column
1603  for (int i = n + 1; i < BYM3; ++i) {
1604  for (int m = 1; m < BXM3; ++m) {
1605  dxytempdy[0][m][i] = dxytempdy[0][m][i + 1];
1606  }
1607  }
1608  }
1609  }
1610 
1611  //combine rows and columns to simulate double pixels
1612 
1613  for (int m = 1; m < BXM3; ++m) {
1614  if (xdouble[m]) {
1615  // Combine the x-rows
1616  for (int n = 1; n < BYM3; ++n) {
1617  dxytempdy[0][m][n] += dxytempdy[0][m + 1][n];
1618  }
1619  // Now shift the remaining pixels over by one row
1620  for (int j = m + 1; j < BXM3; ++j) {
1621  for (int n = 1; n < BYM3; ++n) {
1622  dxytempdy[0][j][n] = dxytempdy[0][j + 1][n];
1623  }
1624  }
1625  }
1626  }
1627 
1628  // Finally, normalize the derivatives and copy the results to the output array
1629 
1630  for (int n = 1; n < BYM3; ++n) {
1631  for (int m = 1; m < BXM3; ++m) {
1632  dpdx2d[1][m][n] = (dxytempdy[1][m][n] - dxytempdy[0][m][n]) / (2. * deltaxy[1]);
1633  }
1634  }
1635  }
1636 
1637  return success_;
1638 } // xytemp
1639 
1640 // *************************************************************************************************************************************
1647 // *************************************************************************************************************************************
1648 
1649 bool SiPixelTemplate2D::xytemp(
1650  float xhit, float yhit, bool ydouble[BYM2], bool xdouble[BXM2], float template2d[BXM2][BYM2]) {
1651  // Interpolate for a new set of track angles
1652 
1653  bool derivatives = false;
1654  float dpdx2d[2][BXM2][BYM2];
1655  float QTemplate;
1656 
1657  return SiPixelTemplate2D::xytemp(xhit, yhit, ydouble, xdouble, template2d, derivatives, dpdx2d, QTemplate);
1658 
1659 } // xytemp
1660 
1661 // *************************************************************************************************************************************
1671 // *************************************************************************************************************************************
1672 
1673 bool SiPixelTemplate2D::xytemp(int id,
1674  float cotalpha,
1675  float cotbeta,
1676  float xhit,
1677  float yhit,
1678  std::vector<bool>& ydouble,
1679  std::vector<bool>& xdouble,
1680  float template2d[BXM2][BYM2]) {
1681  // Local variables
1682 
1683  bool derivatives = false;
1684  float dpdx2d[2][BXM2][BYM2];
1685  float QTemplate;
1686  float locBx = 1.f;
1687  if (cotbeta < 0.f) {
1688  locBx = -1.f;
1689  }
1690  float locBz = locBx;
1691  if (cotalpha < 0.f) {
1692  locBz = -locBx;
1693  }
1694 
1695  bool yd[BYM2], xd[BXM2];
1696 
1697  yd[0] = false;
1698  yd[BYM2 - 1] = false;
1699  for (int i = 0; i < TYSIZE; ++i) {
1700  yd[i + 1] = ydouble[i];
1701  }
1702  xd[0] = false;
1703  xd[BXM2 - 1] = false;
1704  for (int j = 0; j < TXSIZE; ++j) {
1705  xd[j + 1] = xdouble[j];
1706  }
1707 
1708  // Interpolate for a new set of track angles
1709 
1710  if (SiPixelTemplate2D::interpolate(id, cotalpha, cotbeta, locBz, locBx)) {
1711  return SiPixelTemplate2D::xytemp(xhit, yhit, yd, xd, template2d, derivatives, dpdx2d, QTemplate);
1712  } else {
1713  return false;
1714  }
1715 
1716 } // xytemp
1717 
1718 // ************************************************************************************************************
1724 // ************************************************************************************************************
1725 void SiPixelTemplate2D::xysigma2(float qpixel, int index, float& xysig2)
1726 
1727 {
1728  // Interpolate using quantities already stored in the private variables
1729 
1730  // Local variables
1731  float sigi, sigi2, sigi3, sigi4, qscale, err2, err00;
1732 
1733  // Make sure that input is OK
1734 
1735 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1736  if (index < 1 || index >= BYM2) {
1737  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::ysigma2 called with index = " << index << std::endl;
1738  }
1739 #else
1740  assert(index > 0 && index < BYM2);
1741 #endif
1742 
1743  // Define the maximum signal to use in the parameterization
1744 
1745  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1746 
1747  if (qpixel < sxymax_) {
1748  sigi = qpixel;
1749  qscale = 1.f;
1750  } else {
1751  sigi = sxymax_;
1752  qscale = qpixel / sxymax_;
1753  }
1754  sigi2 = sigi * sigi;
1755  sigi3 = sigi2 * sigi;
1756  sigi4 = sigi3 * sigi;
1757  if (index <= T2HYP1) {
1758  err00 = xypary0x0_[0][0] + xypary0x0_[0][1] * sigi + xypary0x0_[0][2] * sigi2 + xypary0x0_[0][3] * sigi3 +
1759  xypary0x0_[0][4] * sigi4;
1760  err2 = err00 +
1761  adcota_ * (xypary0x1_[0][0] + xypary0x1_[0][1] * sigi + xypary0x1_[0][2] * sigi2 + xypary0x1_[0][3] * sigi3 +
1762  xypary0x1_[0][4] * sigi4 - err00) +
1763  adcotb_ * (xypary1x0_[0][0] + xypary1x0_[0][1] * sigi + xypary1x0_[0][2] * sigi2 + xypary1x0_[0][3] * sigi3 +
1764  xypary1x0_[0][4] * sigi4 - err00);
1765  } else {
1766  err00 = xypary0x0_[1][0] + xypary0x0_[1][1] * sigi + xypary0x0_[1][2] * sigi2 + xypary0x0_[1][3] * sigi3 +
1767  xypary0x0_[1][4] * sigi4;
1768  err2 = err00 +
1769  adcota_ * (xypary0x1_[1][0] + xypary0x1_[1][1] * sigi + xypary0x1_[1][2] * sigi2 + xypary0x1_[1][3] * sigi3 +
1770  xypary0x1_[1][4] * sigi4 - err00) +
1771  adcotb_ * (xypary1x0_[1][0] + xypary1x0_[1][1] * sigi + xypary1x0_[1][2] * sigi2 + xypary1x0_[1][3] * sigi3 +
1772  xypary1x0_[1][4] * sigi4 - err00);
1773  }
1774  xysig2 = qscale * err2;
1775  if (xysig2 <= 0.f) {
1776  xysig2 = s50_ * s50_;
1777  }
1778 
1779  return;
1780 
1781 } // End xysigma2
1782 
1783 // ************************************************************************************************************
1785 // ************************************************************************************************************
1786 void SiPixelTemplate2D::landau_par(float lanpar[2][5])
1787 
1788 {
1789  // Interpolate using quantities already stored in the private variables
1790 
1791  for (int i = 0; i < 2; ++i) {
1792  for (int j = 0; j < 5; ++j) {
1793  lanpar[i][j] = lanpar_[i][j];
1794  }
1795  }
1796  return;
1797 
1798 } // End lan_par
SiPixelTemplateEntry2D::pixmax
float pixmax
maximum charge for individual pixels in cluster
Definition: SiPixelTemplate2D.h:48
SiPixelTemplateHeader2D::qscale
float qscale
Charge scaling to match cmssw and pixelav.
Definition: SiPixelTemplate2D.h:83
T2HYP1
#define T2HYP1
Definition: SiPixelTemplateDefs.h:41
mps_fire.i
i
Definition: mps_fire.py:341
SiPixelTemplateEntry2D::qavg
float qavg
average cluster charge for this set of track angles
Definition: SiPixelTemplate2D.h:47
SiPixelTemplateHeader2D::Vbias
float Vbias
detector bias potential in Volts
Definition: SiPixelTemplate2D.h:88
SiPixelTemplateHeader2D::lorxwidth
float lorxwidth
estimate of x-lorentz width for optimal resolution
Definition: SiPixelTemplate2D.h:85
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:1017
SiPixelTemplateEntry2D::runnum
int runnum
< Basic template entry corresponding to a single set of track angles
Definition: SiPixelTemplate2D.h:43
SiPixelTemplateHeader2D::ysize
float ysize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate2D.h:98
SiPixelTemplateStore2D::head
SiPixelTemplateHeader2D head
< template storage structure
Definition: SiPixelTemplate2D.h:103
SiPixelTemplate2D::xysigma2
void xysigma2(float qpixel, int index, float &xysig2)
Definition: SiPixelTemplate2D.cc:1725
TXSIZE
#define TXSIZE
Definition: SiPixelTemplateDefs.h:30
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
T2HY
#define T2HY
Definition: SiPixelTemplateDefs.h:40
SiPixelTemplateEntry2D::clsleny
float clsleny
cluster y-length in pixels at signal height symax/2
Definition: SiPixelTemplate2D.h:61
SiPixelTemplateEntry2D::sxymax
float sxymax
average pixel signal for use of the error parameterization
Definition: SiPixelTemplate2D.h:49
BXM3
#define BXM3
Definition: SiPixelTemplateDefs.h:37
T2HX
#define T2HX
Definition: SiPixelTemplateDefs.h:42
SiPixelTemplateEntry2D::xypar
float xypar[2][5]
pixel uncertainty parameterization
Definition: SiPixelTemplate2D.h:54
SiPixel2DTemplateDBObject::sVector
std::vector< float > sVector() const
Definition: SiPixel2DTemplateDBObject.h:64
ENDL
#define ENDL
Definition: SiPixelTemplate2D.cc:45
SiPixelTemplateEntry2D::scaleyavg
float scaleyavg
average y-error scale factor
Definition: SiPixelTemplate2D.h:67
ysize
const Int_t ysize
Definition: trackSplitPlot.h:43
SiPixelTemplateEntry2D::delyavg
float delyavg
average length difference between template and cluster
Definition: SiPixelTemplate2D.h:68
SiPixelTemplateEntry2D::delysig
float delysig
rms of length difference between template and cluster
Definition: SiPixelTemplate2D.h:69
SiPixelTemplateStore2D::entry
SiPixelTemplateEntry2D ** entry
use 2d entry to store BPix and FPix entries [dynamically allocated]
Definition: SiPixelTemplate2D.h:104
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:31
DeadROC_duringRun.dir
dir
Definition: DeadROC_duringRun.py:23
SiPixelTemplateEntry2D::cotbeta
float cotbeta
cot(beta) is proportional to cluster length in y and is basis of interpolation
Definition: SiPixelTemplate2D.h:45
BXM2
#define BXM2
Definition: SiPixelTemplateDefs.h:36
SiPixelTemplateHeader2D::lorxbias
float lorxbias
estimate of x-lorentz bias
Definition: SiPixelTemplate2D.h:87
SiPixelTemplateEntry2D::iymax
int iymax
the maximum nonzero pixel yindex in template (saves time during interpolation)
Definition: SiPixelTemplate2D.h:51
SiPixelTemplateHeader2D::fbin
float fbin[3]
The QBin definitions in Q_clus/Q_avg.
Definition: SiPixelTemplate2D.h:96
T2XSIZE
#define T2XSIZE
Definition: SiPixelTemplateDefs.h:39
mps_splice.entry
entry
Definition: mps_splice.py:68
SiPixelTemplateHeader2D::ss50
float ss50
1/2 of the single hit dcol threshold in electrons
Definition: SiPixelTemplate2D.h:92
SiPixelTemplateEntry2D::iymin
int iymin
the minimum nonzero pixel yindex in template (saves time during interpolation)
Definition: SiPixelTemplate2D.h:50
SiPixelTemplateEntry2D::xytemp
short int xytemp[7][7][21][7]
templates for y-reconstruction (binned over 1 central pixel)
Definition: SiPixelTemplate2D.h:56
SiPixelTemplateHeader2D::ID
int ID
< template header structure
Definition: SiPixelTemplate2D.h:78
LOGERROR
#define LOGERROR(x)
Definition: SiPixelTemplate2D.cc:43
triggerObjects_cff.id
id
Definition: triggerObjects_cff.py:31
SiPixelTemplateHeader2D::NTxx
int NTxx
number of Template x-entries in each slice
Definition: SiPixelTemplate2D.h:81
SiPixelTemplateEntry2D::mpvvav
float mpvvav
most probable charge in Vavilov distribution (not actually for larger kappa)
Definition: SiPixelTemplate2D.h:63
LOGINFO
#define LOGINFO(x)
Definition: SiPixelTemplate2D.cc:44
visualization-live-secondInstance_cfg.m
m
Definition: visualization-live-secondInstance_cfg.py:60
FrontierConditions_GlobalTag_cff.file
file
Definition: FrontierConditions_GlobalTag_cff.py:13
SiPixelTemplateEntry2D::scalex
float scalex[4]
x-error scale factor in 4 charge bins
Definition: SiPixelTemplate2D.h:70
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
SiPixelTemplateHeader2D::xsize
float xsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate2D.h:97
SiPixelTemplateEntry2D::costrk
float costrk[3]
direction cosines of tracks used to generate this entry
Definition: SiPixelTemplate2D.h:46
SiPixelTemplateEntry2D::offsetx
float offsetx[4]
x-offset in 4 charge bins
Definition: SiPixelTemplate2D.h:72
BYM2
#define BYM2
Definition: SiPixelTemplateDefs.h:28
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
SiPixelTemplateEntry2D::jxmin
int jxmin
the minimum nonzero pixel xindex in template (saves time during interpolation)
Definition: SiPixelTemplate2D.h:52
min
T min(T a, T b)
Definition: MathUtil.h:58
BYM3
#define BYM3
Definition: SiPixelTemplateDefs.h:29
SiPixelTemplateEntry2D::sigmavav
float sigmavav
"sigma" scale fctor for Vavilov distribution
Definition: SiPixelTemplate2D.h:64
SiPixelTemplateEntry2D::cotalpha
float cotalpha
cot(alpha) is proportional to cluster length in x and is basis of interpolation
Definition: SiPixelTemplate2D.h:44
SiPixelTemplate2D::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
Definition: SiPixelTemplate2D.cc:659
SiPixelTemplateHeader2D::title
char title[80]
template title
Definition: SiPixelTemplate2D.h:93
SiPixelTemplateHeader2D::zsize
float zsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate2D.h:99
SiPixelTemplateEntry2D::offsety
float offsety[4]
y-offset in 4 charge bins
Definition: SiPixelTemplate2D.h:73
TYSIZE
#define TYSIZE
Definition: SiPixelTemplateDefs.h:21
SiPixelTemplateEntry2D::chi2scale
float chi2scale
scale factor for the chi2 distribution
Definition: SiPixelTemplate2D.h:58
SiPixelTemplateHeader2D::lorywidth
float lorywidth
estimate of y-lorentz width for optimal resolution
Definition: SiPixelTemplate2D.h:84
SiPixelTemplateEntry2D::clslenx
float clslenx
cluster x-length in pixels at signal height sxmax/2
Definition: SiPixelTemplate2D.h:62
SiPixelTemplateHeader2D::NTy
int NTy
number of Template y entries
Definition: SiPixelTemplate2D.h:79
SiPixelTemplate2D::pushfile
static bool pushfile(int filenum, std::vector< SiPixelTemplateStore2D > &pixelTemp, std::string dir="")
Definition: SiPixelTemplate2D.cc:61
SiPixelTemplateEntry2D::scalexavg
float scalexavg
average x-error scale factor
Definition: SiPixelTemplate2D.h:66
SiPixelTemplateHeader2D::Dtype
int Dtype
detector type (0=BPix, 1=FPix)
Definition: SiPixelTemplate2D.h:82
SiPixelTemplateEntry2D::scaley
float scaley[4]
y-error scale factor in 4 charge bins
Definition: SiPixelTemplate2D.h:71
edm
HLT enums.
Definition: AlignableModifier.h:19
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:1786
edm::FileInPath::fullPath
std::string fullPath() const
Definition: FileInPath.cc:163
SiPixelTemplateEntry2D::chi2ppix
float chi2ppix
average chi^2 per pixel
Definition: SiPixelTemplate2D.h:57
SiPixelTemplateHeader2D::templ_version
int templ_version
Version number of the template to ensure code compatibility.
Definition: SiPixelTemplate2D.h:94
gather_cfg.cout
cout
Definition: gather_cfg.py:144
SiPixelTemplateHeader2D::s50
float s50
1/2 of the multihit dcol threshold in electrons
Definition: SiPixelTemplate2D.h:91
SiPixelTemplateHeader2D::temperature
float temperature
detector temperature in deg K
Definition: SiPixelTemplate2D.h:89
SiPixelTemplateHeader2D::Bfield
float Bfield
Bfield in Tesla.
Definition: SiPixelTemplate2D.h:95
SiPixelTemplate2D::sideload
void sideload(SiPixelTemplateEntry2D *entry, int iDtype, float locBx, float locBz, float lorwdy, float lorwdx, float q50, float fbin[3], float xsize, float ysize, float zsize)
Definition: SiPixelTemplate2D.cc:886
T2YSIZE
#define T2YSIZE
Definition: SiPixelTemplateDefs.h:38
xsize
const Int_t xsize
Definition: trackSplitPlot.h:42
SiPixelTemplateHeader2D::fluence
float fluence
radiation fluence in n_eq/cm^2
Definition: SiPixelTemplate2D.h:90
SiPixelTemplateEntry2D::kappavav
float kappavav
kappa parameter for Vavilov distribution
Definition: SiPixelTemplate2D.h:65
SiPixelTemplateHeader2D::lorybias
float lorybias
estimate of y-lorentz bias
Definition: SiPixelTemplate2D.h:86
SiPixelTemplateHeader2D::NTyx
int NTyx
number of Template y-slices of x entries
Definition: SiPixelTemplate2D.h:80
Generated for CMSSW Reference Manual by doxygen 1.8.11