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