CMS 3D CMS Logo
SiPixelTemplate.cc
Go to the documentation of this file.
1 //
2 // SiPixelTemplate.cc Version 10.21
3 //
4 // Add goodness-of-fit info and spare entries to templates, version number in template header, more error checking
5 // Add correction for (Q_F-Q_L)/(Q_F+Q_L) bias
6 // Add cot(beta) reflection to reduce y-entries and more sophisticated x-interpolation
7 // Fix small index searching bug in interpolate method
8 // Change interpolation indexing to avoid complier complaining about possible un-initialized variables
9 // Replace containers with static arrays in calls to ysigma2 and xsigma2
10 // Add external threshold to calls to ysigma2 and xsigma2, fix parameter signal max for xsigma2
11 // Return to 5 pixel spanning but adjust boundaries to use only when needed
12 // Implement improved (faster) chi2min search that depends on pixel types
13 // Fill template arrays in single calls to this object
14 // Add qmin to the template
15 // Add qscale to match charge scales
16 // Small improvement to x-chisquare interpolation
17 // Enlarge SiPixelTemplateStore to accommodate larger templates and increased alpha acceptance (reduce PT threshold to ~200 MeV)
18 // Change output streams to conform to CMSSW info and error logging.
19 // Store x and y cluster sizes in fractional pixels to facilitate cluster splitting
20 // Add methods to return 3-d templates needed for cluster splitting
21 // Keep interpolated central 9 template bins in private storage and expand/shift in the getter functions (faster for speed=2/3) and easier to build 3d templates
22 // Store error and bias information for the simple chi^2 min position analysis (no interpolation or Q_{FB} corrections) to use in cluster splitting
23 // To save time, the gaussian centers and sigma are not interpolated right now (they aren't currently used). They can be restored by un-commenting lines in the interpolate method.
24 // Add a new method to calculate qbin for input cotbeta and cluster charge. To be used for error estimation of merged clusters in PixelCPEGeneric.
25 // Improve the charge estimation for larger cot(alpha) tracks
26 // Change interpolate method to return false boolean if track angles are outside of range
27 // Add template info and method for truncation information
28 // Change to allow template sizes to be changed at compile time
29 // Fix bug in track angle checking
30 // Accommodate Dave's new DB pushfile which overloads the old method (file input)
31 // Add CPEGeneric error information and expand qbin method to access useful info for PixelCPEGeneric
32 // Fix large cot(alpha) bug in qmin interpolation
33 // Add second qmin to allow a qbin=5 state
34 // Use interpolated chi^2 info for one-pixel clusters
35 // Fix DB pushfile version number checking bug.
36 // Remove assert from qbin method
37 // Replace asserts with exceptions in CMSSW
38 // Change calling sequence to interpolate method to handle cot(beta)<0 for FPix cosmics
39 // Add getter for pixelav Lorentz width estimates to qbin method
40 // Add check on template size to interpolate and qbin methods
41 // Add qbin population information, charge distribution information
42 //
43 //
44 // V7.00 - Decouple BPix and FPix information into separate templates
45 // Add methods to facilitate improved cluster splitting
46 // Fix small charge scaling bug (affects FPix only)
47 // Change y-slice used for the x-template to be closer to the actual cotalpha-cotbeta point
48 // (there is some weak breakdown of x-y factorization in the FPix after irradiation)
49 //
50 //
51 // V8.00 - Add method to calculate a simple 2D template
52 // Reorganize the interpolate method to extract header info only once per ID
53 // V8.01 - Improve simple template normalization
54 // V8.05 - Change qbin normalization to work better after irradiation
55 // V8.10 - Add Vavilov distribution interpolation
56 // V8.11 - Renormalize the x-templates for Guofan's cluster size calculation
57 // V8.12 - Technical fix to qavg issue.
58 // V8.13 - Fix qbin and fastsim interpolaters to avoid changing class variables
59 // V8.20 - Add methods to identify the central pixels in the x- and y-templates (to help align templates with clusters in radiation damaged detectors)
60 // Rename class variables from pxxxx (private xxxx) to xxxx_ to follow standard convention.
61 // Add compiler option to store the template entries in BOOST multiarrays of structs instead of simple c arrays
62 // (allows dynamic resizing template storage and has bounds checking but costs ~10% in execution time).
63 // V8.21 - Add new qbin method to use in cluster splitting
64 // V8.23 - Replace chi-min position errors with merged cluster chi2 probability info
65 // V8.25 - Incorporate VI's speed changes into the current version
66 // V8.26 - Modify the Vavilov lookups to work better with the FPix (offset y-templates)
67 // V8.30 - Change the splitting template generation and access to improve speed and eliminate triple index boost::multiarray
68 // V8.31 - Add correction factor: measured/true charge
69 // V8.31 - Fix version number bug in db object I/O (pushfile)
70 // V8.32 - Check for illegal qmin during loading
71 // V8.33 - Fix small type conversion warnings
72 // V8.40 - Incorporate V.I. optimizations
73 // V9.00 - Expand header to include multi and single dcol thresholds, LA biases, and (variable) Qbin definitions
74 // V9.01 - Protect against negative error squared
75 // V10.00 - Update to work with Phase 1 FPix. Fix some code problems introduced by other maintainers.
76 // V10.01 - Fix initialization style as suggested by S. Krutelyov
77 // V10.10 - Add class variables and methods to be used to correctly calculate the probabilities of single pixel clusters
78 // V10.11 - Allow subdetector ID=5 for FPix R2P2 [allows better internal labeling of templates]
79 // V10.12 - Enforce minimum signal size in pixel charge uncertainty calculation
80 // V10.13 - Update the variable size [SI_PIXEL_TEMPLATE_USE_BOOST] option so that it works with VI's enhancements
81 // V10.20 - Add directory path selection to the ascii pushfile method
82 // V10.21 - Address runtime issues in pushfile() for gcc 7.X due to using tempfile as char string + misc. cleanup [Petar]
83 // V10.22 - Move templateStore to the heap, fix variable name in pushfile() [Petar]
84 
85 
86 // Created by Morris Swartz on 10/27/06.
87 //
88 //
89 
90 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
91 #include <cmath>
92 #else
93 #include <math.h>
94 #endif
95 #include <algorithm>
96 #include <vector>
97 #include "boost/multi_array.hpp"
98 #include <iostream>
99 #include <iomanip>
100 #include <sstream>
101 #include <fstream>
102 #include <list>
103 
104 
105 
106 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
107 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplate.h"
108 #include "CondFormats/SiPixelTransient/interface/SimplePixel.h"
109 #include "FWCore/ParameterSet/interface/FileInPath.h"
110 #include "FWCore/MessageLogger/interface/MessageLogger.h"
111 #include "FWCore/Utilities/interface/isFinite.h"
112 #define LOGERROR(x) LogError(x)
113 #define LOGINFO(x) LogInfo(x)
114 #define LOGWARNING(x) LogWarning(x)
115 #define ENDL " "
116 #include "FWCore/Utilities/interface/Exception.h"
117 using namespace edm;
118 #else
119 #include "SiPixelTemplate.h"
120 #include "SimplePixel.h"
121 #define LOGERROR(x) std::cout << x << ": "
122 #define LOGINFO(x) std::cout << x << ": "
123 #define ENDL std::endl
124 #endif
125 
126 //****************************************************************
131 //****************************************************************
132 bool SiPixelTemplate::pushfile(int filenum, std::vector< SiPixelTemplateStore > & pixelTemp , std::string dir)
133 {
134  // Add template stored in external file numbered filenum to theTemplateStore
135 
136  // Local variables
137  int i, j, k, l;
138  float qavg_avg;
139  char c;
140  const int code_version={17};
141 
142  // Create a filename for this run
143  std::string tempfile = std::to_string(filenum);
144 
145 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
146  // If integer filenum has less than 4 digits, prepend 0's until we have four numerical characters, e.g. "0292"
147  int nzeros = 4-tempfile.length();
148  if (nzeros > 0)
149  tempfile = std::string(nzeros, '0') + tempfile;
151 
152  tempfile = dir + "template_summary_zp" + tempfile + ".out";
153  edm::FileInPath file( tempfile ); // Find the file in CMSSW
154  tempfile = file.fullPath(); // Put it back with the whole path.
155 
156 #else
157  // This is the same as above, but more elegant. (Elegance not allowed in CMSSW...)
158  std::ostringstream tout;
159  tout << "template_summary_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
160  tempfile = tout.str();
161 
162 #endif
163 
164  // Open the template file
165  //
166  std::ifstream in_file(tempfile);
167  if(in_file.is_open() && in_file.good()) {
168 
169 
170  // Create a local template storage entry
171 
172  SiPixelTemplateStore theCurrentTemp;
173 
174  // Read-in a header string first and print it
175 
176  for (i=0; (c=in_file.get()) != '\n'; ++i) {
177  if(i < 79) {theCurrentTemp.head.title[i] = c;}
178  }
179  if(i > 78) {i=78;}
180  theCurrentTemp.head.title[i+1] ='\0';
181  LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
182 
183  // next, the header information
184 
185  in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
186  >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
187  >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
188 
189  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 0A, no template load" << ENDL; return false;}
190 
191  if(theCurrentTemp.head.templ_version > 17) {
192  in_file >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias
193  >> theCurrentTemp.head.lorxbias >> theCurrentTemp.head.fbin[0]
194  >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
195 
196  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 0B, no template load"
197  << ENDL; return false;}
198  } else {
199  theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
200  theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth/2.f;
201  theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth/2.f;
202  theCurrentTemp.head.fbin[0] = 1.5f;
203  theCurrentTemp.head.fbin[1] = 1.00f;
204  theCurrentTemp.head.fbin[2] = 0.85f;
205  }
206 
207  LOGINFO("SiPixelTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
208  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
209  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
210  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
211  << ", 1/2 multi dcol threshold " << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold " << theCurrentTemp.head.ss50
212  << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", y Lorentz Bias " << theCurrentTemp.head.lorybias
213  << ", x Lorentz width " << theCurrentTemp.head.lorxwidth << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
214  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", " << theCurrentTemp.head.fbin[1]
215  << ", " << theCurrentTemp.head.fbin[2]
216  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
217 
218 
219  if(theCurrentTemp.head.templ_version < code_version) {LOGERROR("SiPixelTemplate") << "code expects version " << code_version << " finds "<< theCurrentTemp.head.templ_version <<", no template load" << ENDL; return false;}
220 
221 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
222 
223  // next, layout the 1-d/2-d structures needed to store template
224 
225  theCurrentTemp.cotbetaY = new float[theCurrentTemp.head.NTy];
226  theCurrentTemp.cotbetaX = new float[theCurrentTemp.head.NTyx];
227  theCurrentTemp.cotalphaX = new float[theCurrentTemp.head.NTxx];
228 
229  theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
230 
231  theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
232 
233 #endif
234 
235  // next, loop over all y-angle entries
236 
237  for (i=0; i < theCurrentTemp.head.NTy; ++i) {
238 
239  in_file >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0]
240  >> theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2];
241 
242  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
243 
244  // Calculate the alpha, beta, and cot(beta) for this entry
245 
246  theCurrentTemp.enty[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
247 
248  theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0]/theCurrentTemp.enty[i].costrk[2];
249 
250  theCurrentTemp.enty[i].beta = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
251 
252  theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1]/theCurrentTemp.enty[i].costrk[2];
253 
254  in_file >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].symax >> theCurrentTemp.enty[i].dyone
255  >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].sxmax >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
256 
257  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
258 
259  in_file >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo
260  >> theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clsleny >> theCurrentTemp.enty[i].clslenx;
261 
262  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 3, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
263 
264 
265  if(theCurrentTemp.enty[i].qmin <= 0.) {LOGERROR("SiPixelTemplate") << "Error in template ID " << theCurrentTemp.head.ID << " qmin = " << theCurrentTemp.enty[i].qmin << ", run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
266 
267 
268  for (j=0; j<2; ++j) {
269 
270  in_file >> theCurrentTemp.enty[i].ypar[j][0] >> theCurrentTemp.enty[i].ypar[j][1]
271  >> theCurrentTemp.enty[i].ypar[j][2] >> theCurrentTemp.enty[i].ypar[j][3] >> theCurrentTemp.enty[i].ypar[j][4];
272 
273  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 4, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
274 
275  }
276 
277  for (j=0; j<9; ++j) {
278 
279  for (k=0; k<TYSIZE; ++k) {in_file >> theCurrentTemp.enty[i].ytemp[j][k];}
280 
281  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 5, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
282  }
283 
284  for (j=0; j<2; ++j) {
285 
286  in_file >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1]
287  >> theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
288 
289  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 6, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
290 
291  }
292 
293  qavg_avg = 0.f;
294  for (j=0; j<9; ++j) {
295 
296  for (k=0; k<TXSIZE; ++k) {in_file >> theCurrentTemp.enty[i].xtemp[j][k]; qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];}
297 
298  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 7, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
299  }
300  theCurrentTemp.enty[i].qavg_avg = qavg_avg/9.;
301 
302  for (j=0; j<4; ++j) {
303 
304  in_file >> theCurrentTemp.enty[i].yavg[j] >> theCurrentTemp.enty[i].yrms[j] >> theCurrentTemp.enty[i].ygx0[j] >> theCurrentTemp.enty[i].ygsig[j];
305 
306  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 8, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
307  }
308 
309  for (j=0; j<4; ++j) {
310 
311  in_file >> theCurrentTemp.enty[i].yflpar[j][0] >> theCurrentTemp.enty[i].yflpar[j][1] >> theCurrentTemp.enty[i].yflpar[j][2]
312  >> theCurrentTemp.enty[i].yflpar[j][3] >> theCurrentTemp.enty[i].yflpar[j][4] >> theCurrentTemp.enty[i].yflpar[j][5];
313 
314  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 9, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
315  }
316 
317  for (j=0; j<4; ++j) {
318 
319  in_file >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >> theCurrentTemp.enty[i].xgsig[j];
320 
321  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 10, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
322  }
323 
324  for (j=0; j<4; ++j) {
325 
326  in_file >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >> theCurrentTemp.enty[i].xflpar[j][2]
327  >> theCurrentTemp.enty[i].xflpar[j][3] >> theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
328 
329  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 11, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
330  }
331 
332  for (j=0; j<4; ++j) {
333 
334  in_file >> theCurrentTemp.enty[i].chi2yavg[j] >> theCurrentTemp.enty[i].chi2ymin[j] >> theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j];
335 
336  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 12, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
337  }
338 
339  for (j=0; j<4; ++j) {
340 
341  in_file >> theCurrentTemp.enty[i].yavgc2m[j] >> theCurrentTemp.enty[i].yrmsc2m[j] >> theCurrentTemp.enty[i].chi2yavgc2m[j] >> theCurrentTemp.enty[i].chi2yminc2m[j];
342 
343  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 13, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
344  }
345 
346  for (j=0; j<4; ++j) {
347 
348  in_file >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >> theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
349 
350  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
351  }
352 
353  for (j=0; j<4; ++j) {
354 
355  in_file >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >> theCurrentTemp.enty[i].ygx0gen[j] >> theCurrentTemp.enty[i].ygsiggen[j];
356 
357  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14a, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
358  }
359 
360  for (j=0; j<4; ++j) {
361 
362  in_file >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >> theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
363 
364  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14b, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
365  }
366 
367  in_file >> theCurrentTemp.enty[i].chi2yavgone >> theCurrentTemp.enty[i].chi2yminone >> theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2
368  >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav >> theCurrentTemp.enty[i].r_qMeas_qTrue >> theCurrentTemp.enty[i].spare[0];
369 
370  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 15, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
371 
372  in_file >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >> theCurrentTemp.enty[i].kappavav2 >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1]
373  >> theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracyone >> theCurrentTemp.enty[i].fracxone >> theCurrentTemp.enty[i].fracytwo >> theCurrentTemp.enty[i].fracxtwo;
374  // theCurrentTemp.enty[i].qbfrac[3] = 1. - theCurrentTemp.enty[i].qbfrac[0] - theCurrentTemp.enty[i].qbfrac[1] - theCurrentTemp.enty[i].qbfrac[2];
375 
376  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 16, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
377 
378  }
379 
380  // next, loop over all barrel x-angle entries
381 
382  for (k=0; k < theCurrentTemp.head.NTyx; ++k) {
383 
384  for (i=0; i < theCurrentTemp.head.NTxx; ++i) {
385 
386  in_file >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0]
387  >> theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2];
388 
389  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 17, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
390 
391  // Calculate the alpha, beta, and cot(beta) for this entry
392 
393  theCurrentTemp.entx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
394 
395  theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0]/theCurrentTemp.entx[k][i].costrk[2];
396 
397  theCurrentTemp.entx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
398 
399  theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1]/theCurrentTemp.entx[k][i].costrk[2];
400 
401  in_file >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >> theCurrentTemp.entx[k][i].symax >> theCurrentTemp.entx[k][i].dyone
402  >> theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
403 
404  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 18, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
405 
406  in_file >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >> theCurrentTemp.entx[k][i].dxtwo
407  >> theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clsleny >> theCurrentTemp.entx[k][i].clslenx;
408  // >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
409 
410  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 19, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
411 
412  for (j=0; j<2; ++j) {
413 
414  in_file >> theCurrentTemp.entx[k][i].ypar[j][0] >> theCurrentTemp.entx[k][i].ypar[j][1]
415  >> theCurrentTemp.entx[k][i].ypar[j][2] >> theCurrentTemp.entx[k][i].ypar[j][3] >> theCurrentTemp.entx[k][i].ypar[j][4];
416 
417  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 20, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
418  }
419 
420  for (j=0; j<9; ++j) {
421 
422  for (l=0; l<TYSIZE; ++l) {in_file >> theCurrentTemp.entx[k][i].ytemp[j][l];}
423 
424  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 21, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
425  }
426 
427  for (j=0; j<2; ++j) {
428 
429  in_file >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1]
430  >> theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >> theCurrentTemp.entx[k][i].xpar[j][4];
431 
432 
433  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 22, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
434  }
435 
436  qavg_avg = 0.f;
437  for (j=0; j<9; ++j) {
438 
439  for (l=0; l<TXSIZE; ++l) {in_file >> theCurrentTemp.entx[k][i].xtemp[j][l]; qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];}
440 
441  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 23, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
442  }
443  theCurrentTemp.entx[k][i].qavg_avg = qavg_avg/9.;
444 
445  for (j=0; j<4; ++j) {
446 
447  in_file >> theCurrentTemp.entx[k][i].yavg[j] >> theCurrentTemp.entx[k][i].yrms[j] >> theCurrentTemp.entx[k][i].ygx0[j] >> theCurrentTemp.entx[k][i].ygsig[j];
448 
449  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 24, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
450  }
451 
452  for (j=0; j<4; ++j) {
453 
454  in_file >> theCurrentTemp.entx[k][i].yflpar[j][0] >> theCurrentTemp.entx[k][i].yflpar[j][1] >> theCurrentTemp.entx[k][i].yflpar[j][2]
455  >> theCurrentTemp.entx[k][i].yflpar[j][3] >> theCurrentTemp.entx[k][i].yflpar[j][4] >> theCurrentTemp.entx[k][i].yflpar[j][5];
456 
457  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 25, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
458  }
459 
460  for (j=0; j<4; ++j) {
461 
462  in_file >> theCurrentTemp.entx[k][i].xavg[j] >> theCurrentTemp.entx[k][i].xrms[j] >> theCurrentTemp.entx[k][i].xgx0[j] >> theCurrentTemp.entx[k][i].xgsig[j];
463 
464  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 26, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
465  }
466 
467  for (j=0; j<4; ++j) {
468 
469  in_file >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >> theCurrentTemp.entx[k][i].xflpar[j][2]
470  >> theCurrentTemp.entx[k][i].xflpar[j][3] >> theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
471 
472  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 27, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
473  }
474 
475  for (j=0; j<4; ++j) {
476 
477  in_file >> theCurrentTemp.entx[k][i].chi2yavg[j] >> theCurrentTemp.entx[k][i].chi2ymin[j] >> theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j];
478 
479  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 28, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
480  }
481 
482  for (j=0; j<4; ++j) {
483 
484  in_file >> theCurrentTemp.entx[k][i].yavgc2m[j] >> theCurrentTemp.entx[k][i].yrmsc2m[j] >> theCurrentTemp.entx[k][i].chi2yavgc2m[j] >> theCurrentTemp.entx[k][i].chi2yminc2m[j];
485 
486  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 29, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
487  }
488 
489  for (j=0; j<4; ++j) {
490 
491  in_file >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >> theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
492 
493  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
494  }
495 
496  for (j=0; j<4; ++j) {
497 
498  in_file >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >> theCurrentTemp.entx[k][i].ygx0gen[j] >> theCurrentTemp.entx[k][i].ygsiggen[j];
499 
500  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30a, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
501  }
502 
503  for (j=0; j<4; ++j) {
504 
505  in_file >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j] >> theCurrentTemp.entx[k][i].xgx0gen[j] >> theCurrentTemp.entx[k][i].xgsiggen[j];
506 
507  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30b, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
508  }
509 
510  in_file >> theCurrentTemp.entx[k][i].chi2yavgone >> theCurrentTemp.entx[k][i].chi2yminone >> theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >> theCurrentTemp.entx[k][i].qmin2
511  >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav >> theCurrentTemp.entx[k][i].r_qMeas_qTrue >> theCurrentTemp.entx[k][i].spare[0];
512 
513  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 31, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
514 
515  in_file >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >> theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].qbfrac[0] >> theCurrentTemp.entx[k][i].qbfrac[1]
516  >> theCurrentTemp.entx[k][i].qbfrac[2] >> theCurrentTemp.entx[k][i].fracyone >> theCurrentTemp.entx[k][i].fracxone >> theCurrentTemp.entx[k][i].fracytwo >> theCurrentTemp.entx[k][i].fracxtwo;
517  // theCurrentTemp.entx[k][i].qbfrac[3] = 1. - theCurrentTemp.entx[k][i].qbfrac[0] - theCurrentTemp.entx[k][i].qbfrac[1] - theCurrentTemp.entx[k][i].qbfrac[2];
518 
519  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 32, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
520 
521  }
522  }
523 
524 
525  in_file.close();
526 
527  // Add this template to the store
528 
529  pixelTemp.push_back(theCurrentTemp);
530 
531  postInit(pixelTemp);
532  return true;
533 
534  } else {
535 
536  // If file didn't open, report this
537 
538  LOGERROR("SiPixelTemplate") << "Error opening File" << tempfile << ENDL;
539  return false;
540 
541  }
542 
543 } // TempInit
544 
545 
546 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
547 
548 
549 //****************************************************************
553 //****************************************************************
554 bool SiPixelTemplate::pushfile(const SiPixelTemplateDBObject& dbobject, std::vector< SiPixelTemplateStore > & pixelTemp)
555 {
556  // Add template stored in external dbobject to theTemplateStore
557 
558  // Local variables
559  int i, j, k, l;
560  float qavg_avg;
561  const int code_version={17};
562 
563  // We must create a new object because dbobject must be a const and our stream must not be
564  auto db(dbobject.reader());
565 
566  // Create a local template storage entry
568  auto tmpPtr = std::make_unique<SiPixelTemplateStore>(); // must be allocated on the heap instead
569  auto & theCurrentTemp = *tmpPtr;
570 
571  // Fill the template storage for each template calibration stored in the db
572  for(int m=0; m<db.numOfTempl(); ++m) {
573 
574  // Read-in a header string first and print it
575 
576  SiPixelTemplateDBObject::char2float temp;
577  for (i=0; i<20; ++i) {
578  temp.f = db.sVector()[db.index()];
579  theCurrentTemp.head.title[4*i] = temp.c[0];
580  theCurrentTemp.head.title[4*i+1] = temp.c[1];
581  theCurrentTemp.head.title[4*i+2] = temp.c[2];
582  theCurrentTemp.head.title[4*i+3] = temp.c[3];
583  db.incrementIndex(1);
584  }
585  theCurrentTemp.head.title[79] = '\0';
586  LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
587 
588  // next, the header information
589 
590  db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
591  >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
592  >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
593 
594  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 0A, no template load" << ENDL; return false;}
595 
596  LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title
597  <<" code version = "<<code_version
598  <<" object version "<<theCurrentTemp.head.templ_version
599  << ENDL;
600 
601  if(theCurrentTemp.head.templ_version > 17) {
602  db >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >> theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
603 
604  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 0B, no template load"
605  << ENDL; return false;}
606  } else {
607  theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
608  theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth/2.f;
609  theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth/2.f;
610  theCurrentTemp.head.fbin[0] = 1.50f;
611  theCurrentTemp.head.fbin[1] = 1.00f;
612  theCurrentTemp.head.fbin[2] = 0.85f;
613  //std::cout<<" set fbin "<< theCurrentTemp.head.fbin[0]<<" "<<theCurrentTemp.head.fbin[1]<<" "
614  // <<theCurrentTemp.head.fbin[2]<<std::endl;
615  }
616 
617  LOGINFO("SiPixelTemplate")
618  << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
619  << theCurrentTemp.head.templ_version << ", Bfield = "
620  << theCurrentTemp.head.Bfield<< ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = "
621  << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = "
622  << theCurrentTemp.head.Dtype<< ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
623  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
624  << ", Q-scaling factor " << theCurrentTemp.head.qscale
625  << ", 1/2 multi dcol threshold " << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold "
626  << theCurrentTemp.head.ss50<< ", y Lorentz Width " << theCurrentTemp.head.lorywidth
627  << ", y Lorentz Bias " << theCurrentTemp.head.lorybias
628  << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
629  << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
630  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0]
631  << ", " << theCurrentTemp.head.fbin[1]
632  << ", " << theCurrentTemp.head.fbin[2]
633  << ", pixel x-size " << theCurrentTemp.head.xsize
634  << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
635 
636  if(theCurrentTemp.head.templ_version < code_version) {
637  LOGINFO("SiPixelTemplate") << "code expects version "<< code_version << " finds "
638  << theCurrentTemp.head.templ_version <<", load anyway " << ENDL;
639  //return false; // dk
640  }
641 
642 
643 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
644 
645  // next, layout the 1-d/2-d structures needed to store template
646  theCurrentTemp.cotbetaY = new float[theCurrentTemp.head.NTy];
647  theCurrentTemp.cotbetaX = new float[theCurrentTemp.head.NTyx];
648  theCurrentTemp.cotalphaX = new float[theCurrentTemp.head.NTxx];
649  theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
650  theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
651 
652 #endif
653 
654  // next, loop over all barrel y-angle entries
655 
656  for (i=0; i < theCurrentTemp.head.NTy; ++i) {
657 
658  db >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0]
659  >> theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2];
660 
661  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
662 
663  // Calculate the alpha, beta, and cot(beta) for this entry
664 
665  theCurrentTemp.enty[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
666 
667  theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0]/theCurrentTemp.enty[i].costrk[2];
668 
669  theCurrentTemp.enty[i].beta = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
670 
671  theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1]/theCurrentTemp.enty[i].costrk[2];
672 
673  db >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].symax >> theCurrentTemp.enty[i].dyone
674  >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].sxmax >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
675 
676  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
677 
678  db >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo
679  >> theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clsleny >> theCurrentTemp.enty[i].clslenx;
680  // >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav;
681 
682  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 3, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
683 
684  if(theCurrentTemp.enty[i].qmin <= 0.) {LOGERROR("SiPixelTemplate") << "Error in template ID " << theCurrentTemp.head.ID << " qmin = " << theCurrentTemp.enty[i].qmin << ", run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
685 
686  for (j=0; j<2; ++j) {
687 
688  db >> theCurrentTemp.enty[i].ypar[j][0] >> theCurrentTemp.enty[i].ypar[j][1]
689  >> theCurrentTemp.enty[i].ypar[j][2] >> theCurrentTemp.enty[i].ypar[j][3] >> theCurrentTemp.enty[i].ypar[j][4];
690 
691  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 4, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
692 
693  }
694 
695  for (j=0; j<9; ++j) {
696 
697  for (k=0; k<TYSIZE; ++k) {db >> theCurrentTemp.enty[i].ytemp[j][k];}
698 
699  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 5, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
700  }
701 
702  for (j=0; j<2; ++j) {
703 
704  db >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1]
705  >> theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
706 
707  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 6, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
708 
709  }
710 
711  qavg_avg = 0.f;
712  for (j=0; j<9; ++j) {
713 
714  for (k=0; k<TXSIZE; ++k) {db >> theCurrentTemp.enty[i].xtemp[j][k]; qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];}
715 
716  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 7, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
717  }
718  theCurrentTemp.enty[i].qavg_avg = qavg_avg/9.;
719 
720  for (j=0; j<4; ++j) {
721 
722  db >> theCurrentTemp.enty[i].yavg[j] >> theCurrentTemp.enty[i].yrms[j] >> theCurrentTemp.enty[i].ygx0[j] >> theCurrentTemp.enty[i].ygsig[j];
723 
724  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 8, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
725  }
726 
727  for (j=0; j<4; ++j) {
728 
729  db >> theCurrentTemp.enty[i].yflpar[j][0] >> theCurrentTemp.enty[i].yflpar[j][1] >> theCurrentTemp.enty[i].yflpar[j][2]
730  >> theCurrentTemp.enty[i].yflpar[j][3] >> theCurrentTemp.enty[i].yflpar[j][4] >> theCurrentTemp.enty[i].yflpar[j][5];
731 
732  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 9, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
733  }
734 
735  for (j=0; j<4; ++j) {
736 
737  db >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >> theCurrentTemp.enty[i].xgsig[j];
738 
739  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 10, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
740  }
741 
742  for (j=0; j<4; ++j) {
743 
744  db >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >> theCurrentTemp.enty[i].xflpar[j][2]
745  >> theCurrentTemp.enty[i].xflpar[j][3] >> theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
746 
747  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 11, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
748  }
749 
750  for (j=0; j<4; ++j) {
751 
752  db >> theCurrentTemp.enty[i].chi2yavg[j] >> theCurrentTemp.enty[i].chi2ymin[j] >> theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j];
753 
754  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 12, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
755  }
756 
757  for (j=0; j<4; ++j) {
758 
759  db >> theCurrentTemp.enty[i].yavgc2m[j] >> theCurrentTemp.enty[i].yrmsc2m[j] >> theCurrentTemp.enty[i].chi2yavgc2m[j] >> theCurrentTemp.enty[i].chi2yminc2m[j];
760 
761  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 13, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
762  }
763 
764  for (j=0; j<4; ++j) {
765 
766  db >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >> theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
767 
768  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
769  }
770 
771  for (j=0; j<4; ++j) {
772 
773  db >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >> theCurrentTemp.enty[i].ygx0gen[j] >> theCurrentTemp.enty[i].ygsiggen[j];
774 
775  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14a, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
776  }
777 
778  for (j=0; j<4; ++j) {
779 
780  db >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >> theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
781 
782  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14b, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
783  }
784 
785 
786  db >> theCurrentTemp.enty[i].chi2yavgone >> theCurrentTemp.enty[i].chi2yminone >> theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2
787  >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav >> theCurrentTemp.enty[i].r_qMeas_qTrue >> theCurrentTemp.enty[i].spare[0];
788 
789  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 15, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
790 
791  db >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >> theCurrentTemp.enty[i].kappavav2 >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1]
792  >> theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracyone >> theCurrentTemp.enty[i].fracxone >> theCurrentTemp.enty[i].fracytwo >> theCurrentTemp.enty[i].fracxtwo;
793  // theCurrentTemp.enty[i].qbfrac[3] = 1. - theCurrentTemp.enty[i].qbfrac[0] - theCurrentTemp.enty[i].qbfrac[1] - theCurrentTemp.enty[i].qbfrac[2];
794 
795  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 16, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
796 
797  }
798 
799  // next, loop over all barrel x-angle entries
800 
801  for (k=0; k < theCurrentTemp.head.NTyx; ++k) {
802 
803  for (i=0; i < theCurrentTemp.head.NTxx; ++i) {
804 
805  db >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0]
806  >> theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2];
807 
808  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 17, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
809 
810  // Calculate the alpha, beta, and cot(beta) for this entry
811 
812  theCurrentTemp.entx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
813 
814  theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0]/theCurrentTemp.entx[k][i].costrk[2];
815 
816  theCurrentTemp.entx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
817 
818  theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1]/theCurrentTemp.entx[k][i].costrk[2];
819 
820  db >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >> theCurrentTemp.entx[k][i].symax >> theCurrentTemp.entx[k][i].dyone
821  >> theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
822 
823  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 18, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
824 
825  db >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >> theCurrentTemp.entx[k][i].dxtwo
826  >> theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clsleny >> theCurrentTemp.entx[k][i].clslenx;
827  // >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
828 
829  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 19, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
830 
831  for (j=0; j<2; ++j) {
832 
833  db >> theCurrentTemp.entx[k][i].ypar[j][0] >> theCurrentTemp.entx[k][i].ypar[j][1]
834  >> theCurrentTemp.entx[k][i].ypar[j][2] >> theCurrentTemp.entx[k][i].ypar[j][3] >> theCurrentTemp.entx[k][i].ypar[j][4];
835 
836  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 20, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
837  }
838 
839  for (j=0; j<9; ++j) {
840 
841  for (l=0; l<TYSIZE; ++l) {db >> theCurrentTemp.entx[k][i].ytemp[j][l];}
842 
843  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 21, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
844  }
845 
846 
847 
848  for (j=0; j<2; ++j) {
849 
850  db >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1]
851  >> theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >> theCurrentTemp.entx[k][i].xpar[j][4];
852 
853 
854  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 22, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
855  }
856 
857  qavg_avg = 0.f;
858  for (j=0; j<9; ++j) {
859 
860  for (l=0; l<TXSIZE; ++l) {db >> theCurrentTemp.entx[k][i].xtemp[j][l]; qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];}
861 
862  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 23, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
863  }
864  theCurrentTemp.entx[k][i].qavg_avg = qavg_avg/9.;
865 
866  for (j=0; j<4; ++j) {
867 
868  db >> theCurrentTemp.entx[k][i].yavg[j] >> theCurrentTemp.entx[k][i].yrms[j] >> theCurrentTemp.entx[k][i].ygx0[j] >> theCurrentTemp.entx[k][i].ygsig[j];
869 
870  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 24, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
871  }
872 
873  for (j=0; j<4; ++j) {
874 
875  db >> theCurrentTemp.entx[k][i].yflpar[j][0] >> theCurrentTemp.entx[k][i].yflpar[j][1] >> theCurrentTemp.entx[k][i].yflpar[j][2]
876  >> theCurrentTemp.entx[k][i].yflpar[j][3] >> theCurrentTemp.entx[k][i].yflpar[j][4] >> theCurrentTemp.entx[k][i].yflpar[j][5];
877 
878  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 25, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
879  }
880 
881  for (j=0; j<4; ++j) {
882 
883  db >> theCurrentTemp.entx[k][i].xavg[j] >> theCurrentTemp.entx[k][i].xrms[j] >> theCurrentTemp.entx[k][i].xgx0[j] >> theCurrentTemp.entx[k][i].xgsig[j];
884 
885  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 26, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
886  }
887 
888  for (j=0; j<4; ++j) {
889 
890  db >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >> theCurrentTemp.entx[k][i].xflpar[j][2]
891  >> theCurrentTemp.entx[k][i].xflpar[j][3] >> theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
892 
893  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 27, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
894  }
895 
896  for (j=0; j<4; ++j) {
897 
898  db >> theCurrentTemp.entx[k][i].chi2yavg[j] >> theCurrentTemp.entx[k][i].chi2ymin[j] >> theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j];
899 
900  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 28, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
901  }
902 
903  for (j=0; j<4; ++j) {
904 
905  db >> theCurrentTemp.entx[k][i].yavgc2m[j] >> theCurrentTemp.entx[k][i].yrmsc2m[j] >> theCurrentTemp.entx[k][i].chi2yavgc2m[j] >> theCurrentTemp.entx[k][i].chi2yminc2m[j];
906 
907  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 29, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
908  }
909 
910  for (j=0; j<4; ++j) {
911 
912  db >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >> theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
913 
914  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
915  }
916 
917  for (j=0; j<4; ++j) {
918 
919  db >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >> theCurrentTemp.entx[k][i].ygx0gen[j] >> theCurrentTemp.entx[k][i].ygsiggen[j];
920 
921  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30a, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
922  }
923 
924  for (j=0; j<4; ++j) {
925 
926  db >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j] >> theCurrentTemp.entx[k][i].xgx0gen[j] >> theCurrentTemp.entx[k][i].xgsiggen[j];
927 
928  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30b, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
929  }
930 
931 
932  db >> theCurrentTemp.entx[k][i].chi2yavgone >> theCurrentTemp.entx[k][i].chi2yminone >> theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >> theCurrentTemp.entx[k][i].qmin2
933  >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav >> theCurrentTemp.entx[k][i].r_qMeas_qTrue >> theCurrentTemp.entx[k][i].spare[0];
934 
935  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 31, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
936 
937  db >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >> theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].qbfrac[0] >> theCurrentTemp.entx[k][i].qbfrac[1]
938  >> theCurrentTemp.entx[k][i].qbfrac[2] >> theCurrentTemp.entx[k][i].fracyone >> theCurrentTemp.entx[k][i].fracxone >> theCurrentTemp.entx[k][i].fracytwo >> theCurrentTemp.entx[k][i].fracxtwo;
939  // theCurrentTemp.entx[k][i].qbfrac[3] = 1. - theCurrentTemp.entx[k][i].qbfrac[0] - theCurrentTemp.entx[k][i].qbfrac[1] - theCurrentTemp.entx[k][i].qbfrac[2];
940 
941  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 32, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
942 
943  }
944  }
945 
946 
947  // Add this template to the store
948 
949  pixelTemp.push_back(theCurrentTemp);
950 
951  }
952  postInit(pixelTemp);
953  return true;
954 
955 } // TempInit
956 
957 #endif
958 
959 
960 void SiPixelTemplate::postInit(std::vector< SiPixelTemplateStore > & thePixelTemp_) {
961  /*
962  std::cout << "SiPixelTemplate size " << thePixelTemp_.size() << std::endl;
963  #ifndef SI_PIXEL_TEMPLATE_USE_BOOST
964  std::cout <<"uses C arrays" << std::endl;
965  #endif
966 
967  int i=0;
968  for (auto & templ : thePixelTemp_) {
969  std::cout << i <<':' << templ.head.ID << ' ' << templ.head.NTy <<','<< templ.head.NTyx <<','<< templ.head.NTxx << std::endl;
970  for ( auto iy=1; iy<templ.head.NTy; ++iy ) { auto & ent = templ.enty[iy]; std::cout << ent.cotbeta <<',' << ent.cotbeta-templ.enty[iy-1].cotbeta << ' '; }
971  std::cout << std::endl;
972  for ( auto ix=1; ix<templ.head.NTxx; ++ix ){ auto & ent = templ.entx[0][ix]; std::cout << ent.cotalpha <<','<< ent.cotalpha-templ.entx[0][ix-1].cotalpha << ' ';}
973  std::cout << std::endl;
974  ++i;
975  }
976  */
977 
978  for (auto & templ : thePixelTemp_) {
979  for ( auto iy=0; iy<templ.head.NTy; ++iy ) templ.cotbetaY[iy]=templ.enty[iy].cotbeta;
980  for ( auto iy=0; iy<templ.head.NTyx; ++iy ) templ.cotbetaX[iy]= templ.entx[iy][0].cotbeta;
981  for ( auto ix=0; ix<templ.head.NTxx; ++ix ) templ.cotalphaX[ix]=templ.entx[0][ix].cotalpha;
982  }
983 
984 }
985 
986 
987 
988 
989 // ************************************************************************************************************
1000 // ************************************************************************************************************
1001 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx){
1002  // Interpolate for a new set of track angles
1003 
1004  // Local variables
1005  int i, j;
1006  int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx;
1007  float yratio, yxratio, xxratio, sxmax, qcorrect, qxtempcor, symax, chi2xavgone, chi2xminone, cota, cotb, cotalpha0, cotbeta0;
1008  bool flip_x, flip_y;
1009  // std::vector <float> xrms(4), xgsig(4), xrmsc2m(4);
1010  float chi2xavg[4], chi2xmin[4], chi2xavgc2m[4], chi2xminc2m[4];
1011 
1012 
1013  // Check to see if interpolation is valid
1014  if(id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
1015 
1016  cota_current_ = cotalpha; cotb_current_ = cotbeta; success_ = true;
1017 
1018  if(id != id_current_) {
1019 
1020  // Find the index corresponding to id
1021 
1022  index_id_ = -1;
1023  for(i=0; i<(int)thePixelTemp_.size(); ++i) {
1024 
1025  //std::cout<<i<<" "<<id<<" "<<thePixelTemp_[i].head.ID<<std::endl;
1026 
1027  if(id == thePixelTemp_[i].head.ID) {
1028 
1029  index_id_ = i;
1030  id_current_ = id;
1031 
1032  // Copy the charge scaling factor to the private variable
1033 
1034  qscale_ = thePixelTemp_[index_id_].head.qscale;
1035 
1036  // Copy the pseudopixel signal size to the private variable
1037 
1038  s50_ = thePixelTemp_[index_id_].head.s50;
1039 
1040  // Copy Qbinning info to private variables
1041 
1042  for(j=0; j<3; ++j) {fbin_[j] = thePixelTemp_[index_id_].head.fbin[j];}
1043  //std::cout<<" set fbin "<< fbin_[0]<<" "<<fbin_[1]<<" "<<fbin_[2]<<std::endl;
1044 
1045  // Pixel sizes to the private variables
1046 
1047  xsize_ = thePixelTemp_[index_id_].head.xsize;
1048  ysize_ = thePixelTemp_[index_id_].head.ysize;
1049  zsize_ = thePixelTemp_[index_id_].head.zsize;
1050 
1051  break;
1052  }
1053  }
1054  }
1055 
1056 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1057  if(index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
1058  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::interpolate can't find needed template ID = " << id << std::endl;
1059  }
1060 
1061  //check for nan's
1062  if (!edm::isFinite(cotalpha) || !edm::isFinite(cotbeta)) {
1063  success_ = false;
1064  return success_;
1065  }
1066 #else
1067  assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
1068 #endif
1069 
1070  // Interpolate the absolute value of cot(beta)
1071 
1072  abs_cotb_ = std::abs(cotbeta);
1073 
1074  // qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
1075 
1076  cotalpha0 = thePixelTemp_[index_id_].enty[0].cotalpha;
1077  qcorrect=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta*cotbeta+cotalpha0*cotalpha0));
1078 
1079  // for some cosmics, the ususal gymnastics are incorrect
1080  cota = cotalpha;
1081  cotb = abs_cotb_;
1082  flip_x = false;
1083  flip_y = false;
1084  switch(thePixelTemp_[index_id_].head.Dtype) {
1085  case 0:
1086  if(cotbeta < 0.f) {flip_y = true;}
1087  break;
1088  case 1:
1089  if(locBz < 0.f) {
1090  cotb = cotbeta;
1091  } else {
1092  cotb = -cotbeta;
1093  flip_y = true;
1094  }
1095  break;
1096  case 2:
1097  case 3:
1098  case 4:
1099  case 5:
1100  if(locBx*locBz < 0.f) {
1101  cota = -cotalpha;
1102  flip_x = true;
1103  }
1104  if(locBx > 0.f) {
1105  cotb = cotbeta;
1106  } else {
1107  cotb = -cotbeta;
1108  flip_y = true;
1109  }
1110  break;
1111  default:
1112 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1113  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
1114 #else
1115  std::cout << "SiPixelTemplate::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
1116 #endif
1117  }
1118 
1119  Ny = thePixelTemp_[index_id_].head.NTy;
1120  Nyx = thePixelTemp_[index_id_].head.NTyx;
1121  Nxx = thePixelTemp_[index_id_].head.NTxx;
1122 
1123 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1124  if(Ny < 2 || Nyx < 1 || Nxx < 2) {
1125  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
1126  }
1127 #else
1128  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
1129 #endif
1130  imaxx = Nyx - 1;
1131  imidy = Nxx/2;
1132 
1133  // next, loop over all y-angle entries
1134 
1135  ilow = 0;
1136  yratio = 0.f;
1137 
1138  if(cotb >= thePixelTemp_[index_id_].enty[Ny-1].cotbeta) {
1139 
1140  ilow = Ny-2;
1141  yratio = 1.;
1142  success_ = false;
1143 
1144  } else {
1145 
1146  if(cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
1147 
1148  for (i=0; i<Ny-1; ++i) {
1149 
1150  if( thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i+1].cotbeta) {
1151 
1152  ilow = i;
1153  yratio = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta)/(thePixelTemp_[index_id_].enty[i+1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
1154  break;
1155  }
1156  }
1157  } else { success_ = false; }
1158  }
1159 
1160  ihigh=ilow + 1;
1161 
1162  // Interpolate/store all y-related quantities (flip displacements when flip_y)
1163 
1164  yratio_ = yratio;
1165  qavg_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].qavg + yratio*thePixelTemp_[index_id_].enty[ihigh].qavg;
1166  qavg_ *= qcorrect;
1167  symax = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].symax + yratio*thePixelTemp_[index_id_].enty[ihigh].symax;
1168  syparmax_ = symax;
1169  sxmax = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sxmax + yratio*thePixelTemp_[index_id_].enty[ihigh].sxmax;
1170  dyone_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].dyone + yratio*thePixelTemp_[index_id_].enty[ihigh].dyone;
1171  if(flip_y) {dyone_ = -dyone_;}
1172  syone_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].syone + yratio*thePixelTemp_[index_id_].enty[ihigh].syone;
1173  dytwo_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].dytwo + yratio*thePixelTemp_[index_id_].enty[ihigh].dytwo;
1174  if(flip_y) {dytwo_ = -dytwo_;}
1175  sytwo_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sytwo + yratio*thePixelTemp_[index_id_].enty[ihigh].sytwo;
1176  qmin_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].qmin + yratio*thePixelTemp_[index_id_].enty[ihigh].qmin;
1177  qmin_ *= qcorrect;
1178  qmin2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].qmin2 + yratio*thePixelTemp_[index_id_].enty[ihigh].qmin2;
1179  qmin2_ *= qcorrect;
1180  mpvvav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].mpvvav + yratio*thePixelTemp_[index_id_].enty[ihigh].mpvvav;
1181  mpvvav_ *= qcorrect;
1182  sigmavav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sigmavav + yratio*thePixelTemp_[index_id_].enty[ihigh].sigmavav;
1183  kappavav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].kappavav + yratio*thePixelTemp_[index_id_].enty[ihigh].kappavav;
1184  mpvvav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].mpvvav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].mpvvav2;
1185  mpvvav2_ *= qcorrect;
1186  sigmavav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sigmavav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].sigmavav2;
1187  kappavav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].kappavav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].kappavav2;
1188  clsleny_ = fminf(thePixelTemp_[index_id_].enty[ilow].clsleny, thePixelTemp_[index_id_].enty[ihigh].clsleny);
1189  qavg_avg_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].qavg_avg + yratio*thePixelTemp_[index_id_].enty[ihigh].qavg_avg;
1190  qavg_avg_ *= qcorrect;
1191  for(i=0; i<2 ; ++i) {
1192  for(j=0; j<5 ; ++j) {
1193  // Charge loss switches sides when cot(beta) changes sign
1194  if(flip_y) {
1195  yparl_[1-i][j] = thePixelTemp_[index_id_].enty[ilow].ypar[i][j];
1196  yparh_[1-i][j] = thePixelTemp_[index_id_].enty[ihigh].ypar[i][j];
1197  } else {
1198  yparl_[i][j] = thePixelTemp_[index_id_].enty[ilow].ypar[i][j];
1199  yparh_[i][j] = thePixelTemp_[index_id_].enty[ihigh].ypar[i][j];
1200  }
1201  if(flip_x) {
1202  xparly0_[1-i][j] = thePixelTemp_[index_id_].enty[ilow].xpar[i][j];
1203  xparhy0_[1-i][j] = thePixelTemp_[index_id_].enty[ihigh].xpar[i][j];
1204  } else {
1205  xparly0_[i][j] = thePixelTemp_[index_id_].enty[ilow].xpar[i][j];
1206  xparhy0_[i][j] = thePixelTemp_[index_id_].enty[ihigh].xpar[i][j];
1207  }
1208  }
1209  }
1210 
1211  for(i=0; i<4; ++i) {
1212  yavg_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].yavg[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].yavg[i];
1213  if(flip_y) {yavg_[i] = -yavg_[i];}
1214  yrms_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].yrms[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].yrms[i];
1215  // ygx0_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].ygx0[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].ygx0[i];
1216  // if(flip_y) {ygx0_[i] = -ygx0_[i];}
1217  // ygsig_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].ygsig[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].ygsig[i];
1218  // xrms[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].xrms[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].xrms[i];
1219  // xgsig[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].xgsig[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].xgsig[i];
1220  chi2yavg_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2yavg[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2yavg[i];
1221  chi2ymin_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2ymin[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2ymin[i];
1222  chi2xavg[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xavg[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xavg[i];
1223  chi2xmin[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xmin[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xmin[i];
1224  yavgc2m_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].yavgc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].yavgc2m[i];
1225  if(flip_y) {yavgc2m_[i] = -yavgc2m_[i];}
1226  yrmsc2m_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].yrmsc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].yrmsc2m[i];
1227  chi2yavgc2m_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2yavgc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2yavgc2m[i];
1228  // if(flip_y) {chi2yavgc2m_[i] = -chi2yavgc2m_[i];}
1229  chi2yminc2m_[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2yminc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2yminc2m[i];
1230  // xrmsc2m[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].xrmsc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].xrmsc2m[i];
1231  chi2xavgc2m[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xavgc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xavgc2m[i];
1232  chi2xminc2m[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xminc2m[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xminc2m[i];
1233  for(j=0; j<6 ; ++j) {
1234  yflparl_[i][j] = thePixelTemp_[index_id_].enty[ilow].yflpar[i][j];
1235  yflparh_[i][j] = thePixelTemp_[index_id_].enty[ihigh].yflpar[i][j];
1236 
1237  // Since Q_fl is odd under cotbeta, it flips qutomatically, change only even terms
1238 
1239  if(flip_y && (j == 0 || j == 2 || j == 4)) {
1240  yflparl_[i][j] = - yflparl_[i][j];
1241  yflparh_[i][j] = - yflparh_[i][j];
1242  }
1243  }
1244  }
1245 
1247 
1248  chi2yavgone_=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2yavgone + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2yavgone;
1249  chi2yminone_=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2yminone + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2yminone;
1250  chi2xavgone=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xavgone + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xavgone;
1251  chi2xminone=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].chi2xminone + yratio*thePixelTemp_[index_id_].enty[ihigh].chi2xminone;
1252 
1253  fracyone_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].fracyone + yratio*thePixelTemp_[index_id_].enty[ihigh].fracyone;
1254  fracytwo_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].fracytwo + yratio*thePixelTemp_[index_id_].enty[ihigh].fracytwo;
1255  // for(i=0; i<10; ++i) {
1256  // pyspare[i]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].yspare[i] + yratio*thePixelTemp_[index_id_].enty[ihigh].yspare[i];
1257  // }
1258 
1259  // Interpolate and build the y-template
1260 
1261  for(i=0; i<9; ++i) {
1262  ytemp_[i][0] = 0.f;
1263  ytemp_[i][1] = 0.f;
1264  ytemp_[i][BYM2] = 0.f;
1265  ytemp_[i][BYM1] = 0.f;
1266  for(j=0; j<TYSIZE; ++j) {
1267 
1268  // Flip the basic y-template when the cotbeta is negative
1269 
1270  if(flip_y) {
1271  ytemp_[8-i][BYM3-j]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].ytemp[i][j] + yratio*thePixelTemp_[index_id_].enty[ihigh].ytemp[i][j];
1272  } else {
1273  ytemp_[i][j+2]=(1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].ytemp[i][j] + yratio*thePixelTemp_[index_id_].enty[ihigh].ytemp[i][j];
1274  }
1275  }
1276  }
1277 
1278  // next, loop over all x-angle entries, first, find relevant y-slices
1279 
1280  iylow = 0;
1281  yxratio = 0.f;
1282 
1283  if(abs_cotb_ >= thePixelTemp_[index_id_].entx[Nyx-1][0].cotbeta) {
1284 
1285  iylow = Nyx-2;
1286  yxratio = 1.f;
1287 
1288  } else if(abs_cotb_ >= thePixelTemp_[index_id_].entx[0][0].cotbeta) {
1289 
1290  for (i=0; i<Nyx-1; ++i) {
1291 
1292  if( thePixelTemp_[index_id_].entx[i][0].cotbeta <= abs_cotb_ && abs_cotb_ < thePixelTemp_[index_id_].entx[i+1][0].cotbeta) {
1293 
1294  iylow = i;
1295  yxratio = (abs_cotb_ - thePixelTemp_[index_id_].entx[i][0].cotbeta)/(thePixelTemp_[index_id_].entx[i+1][0].cotbeta - thePixelTemp_[index_id_].entx[i][0].cotbeta);
1296  break;
1297  }
1298  }
1299  }
1300 
1301  iyhigh=iylow + 1;
1302 
1303  ilow = 0;
1304  xxratio = 0.f;
1305 
1306  if(cota >= thePixelTemp_[index_id_].entx[0][Nxx-1].cotalpha) {
1307 
1308  ilow = Nxx-2;
1309  xxratio = 1.f;
1310  success_ = false;
1311 
1312  } else {
1313 
1314  if(cota >= thePixelTemp_[index_id_].entx[0][0].cotalpha) {
1315 
1316  for (i=0; i<Nxx-1; ++i) {
1317 
1318  if( thePixelTemp_[index_id_].entx[0][i].cotalpha <= cota && cota < thePixelTemp_[index_id_].entx[0][i+1].cotalpha) {
1319 
1320  ilow = i;
1321  xxratio = (cota - thePixelTemp_[index_id_].entx[0][i].cotalpha)/(thePixelTemp_[index_id_].entx[0][i+1].cotalpha - thePixelTemp_[index_id_].entx[0][i].cotalpha);
1322  break;
1323  }
1324  }
1325  } else { success_ = false; }
1326  }
1327 
1328  ihigh=ilow + 1;
1329 
1330  // Interpolate/store all x-related quantities
1331 
1332  yxratio_ = yxratio;
1333  xxratio_ = xxratio;
1334 
1335  // sxparmax defines the maximum charge for which the parameters xpar are defined (not rescaled by cotbeta)
1336 
1337  sxparmax_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].sxmax + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].sxmax;
1338  sxmax_ = sxparmax_;
1339  if(thePixelTemp_[index_id_].entx[imaxx][imidy].sxmax != 0.f) {sxmax_=sxmax_/thePixelTemp_[index_id_].entx[imaxx][imidy].sxmax*sxmax;}
1340  symax_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].symax + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].symax;
1341  if(thePixelTemp_[index_id_].entx[imaxx][imidy].symax != 0.f) {symax_=symax_/thePixelTemp_[index_id_].entx[imaxx][imidy].symax*symax;}
1342  dxone_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[0][ilow].dxone + xxratio*thePixelTemp_[index_id_].entx[0][ihigh].dxone;
1343  if(flip_x) {dxone_ = -dxone_;}
1344  sxone_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[0][ilow].sxone + xxratio*thePixelTemp_[index_id_].entx[0][ihigh].sxone;
1345  dxtwo_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[0][ilow].dxtwo + xxratio*thePixelTemp_[index_id_].entx[0][ihigh].dxtwo;
1346  if(flip_x) {dxtwo_ = -dxtwo_;}
1347  sxtwo_ = (1.f - xxratio)*thePixelTemp_[index_id_].entx[0][ilow].sxtwo + xxratio*thePixelTemp_[index_id_].entx[0][ihigh].sxtwo;
1348  clslenx_ = fminf(thePixelTemp_[index_id_].entx[0][ilow].clslenx, thePixelTemp_[index_id_].entx[0][ihigh].clslenx);
1349 
1350  for(i=0; i<2 ; ++i) {
1351  for(j=0; j<5 ; ++j) {
1352  // Charge loss switches sides when cot(alpha) changes sign
1353  if(flip_x) {
1354  xpar0_[1-i][j] = thePixelTemp_[index_id_].entx[imaxx][imidy].xpar[i][j];
1355  xparl_[1-i][j] = thePixelTemp_[index_id_].entx[imaxx][ilow].xpar[i][j];
1356  xparh_[1-i][j] = thePixelTemp_[index_id_].entx[imaxx][ihigh].xpar[i][j];
1357  } else {
1358  xpar0_[i][j] = thePixelTemp_[index_id_].entx[imaxx][imidy].xpar[i][j];
1359  xparl_[i][j] = thePixelTemp_[index_id_].entx[imaxx][ilow].xpar[i][j];
1360  xparh_[i][j] = thePixelTemp_[index_id_].entx[imaxx][ihigh].xpar[i][j];
1361  }
1362  }
1363  }
1364 
1365  // pixmax is the maximum allowed pixel charge (used for truncation)
1366 
1367  pixmax_=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].pixmax + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].pixmax)
1368  +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].pixmax + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].pixmax);
1369 
1370 
1371  r_qMeas_qTrue_=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].r_qMeas_qTrue + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].r_qMeas_qTrue)
1372  +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].r_qMeas_qTrue + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].r_qMeas_qTrue);
1373 
1374 
1375  for(i=0; i<4; ++i) {
1376  xavg_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xavg[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xavg[i])
1377  +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xavg[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xavg[i]);
1378  if(flip_x) {xavg_[i] = -xavg_[i];}
1379 
1380  xrms_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xrms[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xrms[i])
1381  +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xrms[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xrms[i]);
1382 
1383  // xgx0_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xgx0[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xgx0[i])
1384  // +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xgx0[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xgx0[i]);
1385 
1386  // xgsig_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xgsig[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xgsig[i])
1387  // +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xgsig[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xgsig[i]);
1388 
1389  xavgc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xavgc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xavgc2m[i])
1390  +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xavgc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xavgc2m[i]);
1391  if(flip_x) {xavgc2m_[i] = -xavgc2m_[i];}
1392 
1393  xrmsc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xrmsc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xrmsc2m[i])
1394  +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xrmsc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xrmsc2m[i]);
1395 
1396  // chi2xavgc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].chi2xavgc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].chi2xavgc2m[i])
1397  // +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xavgc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xavgc2m[i]);
1398 
1399  // chi2xminc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].chi2xminc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].chi2xminc2m[i])
1400  // +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xminc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xminc2m[i]);
1401  //
1402  // Try new interpolation scheme
1403  //
1404  // chi2xavg_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].chi2xavg[i] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].chi2xavg[i]);
1405  // if(thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i] != 0.f) {chi2xavg_[i]=chi2xavg_[i]/thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i]*chi2xavg[i];}
1406  //
1407  // chi2xmin_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].chi2xmin[i] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].chi2xmin[i]);
1408  // if(thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i] != 0.f) {chi2xmin_[i]=chi2xmin_[i]/thePixelTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i]*chi2xmin[i];}
1409  //
1410  chi2xavg_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xavg[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xavg[i]);
1411  if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavg[i] != 0.f) {chi2xavg_[i]=chi2xavg_[i]/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavg[i]*chi2xavg[i];}
1412 
1413  chi2xmin_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xmin[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xmin[i]);
1414  if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xmin[i] != 0.f) {chi2xmin_[i]=chi2xmin_[i]/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xmin[i]*chi2xmin[i];}
1415 
1416  chi2xavgc2m_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xavgc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xavgc2m[i]);
1417  if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavgc2m[i] != 0.f) {chi2xavgc2m_[i]=chi2xavgc2m_[i]/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavgc2m[i]*chi2xavgc2m[i];}
1418 
1419  chi2xminc2m_[i]=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xminc2m[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xminc2m[i]);
1420  if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xminc2m[i] != 0.f) {chi2xminc2m_[i]=chi2xminc2m_[i]/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xminc2m[i]*chi2xminc2m[i];}
1421 
1422  for(j=0; j<6 ; ++j) {
1423  xflparll_[i][j] = thePixelTemp_[index_id_].entx[iylow][ilow].xflpar[i][j];
1424  xflparlh_[i][j] = thePixelTemp_[index_id_].entx[iylow][ihigh].xflpar[i][j];
1425  xflparhl_[i][j] = thePixelTemp_[index_id_].entx[iyhigh][ilow].xflpar[i][j];
1426  xflparhh_[i][j] = thePixelTemp_[index_id_].entx[iyhigh][ihigh].xflpar[i][j];
1427  // Since Q_fl is odd under cotalpha, it flips qutomatically, change only even terms
1428  if(flip_x && (j == 0 || j == 2 || j == 4)) {
1429  xflparll_[i][j] = -xflparll_[i][j];
1430  xflparlh_[i][j] = -xflparlh_[i][j];
1431  xflparhl_[i][j] = -xflparhl_[i][j];
1432  xflparhh_[i][j] = -xflparhh_[i][j];
1433  }
1434  }
1435  }
1436 
1437 
1438  // Do the spares next
1439 
1440  chi2xavgone_=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xavgone + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xavgone);
1441  if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavgone != 0.f) {chi2xavgone_=chi2xavgone_/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xavgone*chi2xavgone;}
1442 
1443  chi2xminone_=((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].chi2xminone + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].chi2xminone);
1444  if(thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xminone != 0.f) {chi2xminone_=chi2xminone_/thePixelTemp_[index_id_].entx[iyhigh][imidy].chi2xminone*chi2xminone;}
1445 
1446  fracxone_ = (1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].fracxone + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].fracxone)
1447  +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].fracxone + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].fracxone);
1448  fracxtwo_ = (1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].fracxtwo + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].fracxtwo)
1449  +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].fracxtwo + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].fracxtwo);
1450 
1451  // for(i=0; i<10; ++i) {
1452  // pxspare[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iylow][ilow].xspare[i] + xxratio*thePixelTemp_[index_id_].entx[iylow][ihigh].xspare[i])
1453  // +yxratio*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xspare[i] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xspare[i]);
1454  // }
1455 
1456  // Interpolate and build the x-template
1457 
1458  // qxtempcor corrects the total charge to the actual track angles (not actually needed for the template fits, but useful for Guofan)
1459 
1460  cotbeta0 = thePixelTemp_[index_id_].entx[iyhigh][0].cotbeta;
1461  qxtempcor=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta0*cotbeta0+cotalpha*cotalpha));
1462 
1463  for(i=0; i<9; ++i) {
1464  xtemp_[i][0] = 0.f;
1465  xtemp_[i][1] = 0.f;
1466  xtemp_[i][BXM2] = 0.f;
1467  xtemp_[i][BXM1] = 0.f;
1468  for(j=0; j<TXSIZE; ++j) {
1469  // Take next largest x-slice for the x-template (it reduces bias in the forward direction after irradiation)
1470  // xtemp_[i][j+2]=(1.f - xxratio)*thePixelTemp_[index_id_].entx[imaxx][ilow].xtemp[i][j] + xxratio*thePixelTemp_[index_id_].entx[imaxx][ihigh].xtemp[i][j];
1471  // xtemp_[i][j+2]=(1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xtemp[i][j] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xtemp[i][j];
1472  if(flip_x) {
1473  xtemp_[8-i][BXM3-j]=qxtempcor*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xtemp[i][j] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xtemp[i][j]);
1474  } else {
1475  xtemp_[i][j+2]=qxtempcor*((1.f - xxratio)*thePixelTemp_[index_id_].entx[iyhigh][ilow].xtemp[i][j] + xxratio*thePixelTemp_[index_id_].entx[iyhigh][ihigh].xtemp[i][j]);
1476  }
1477  }
1478  }
1479 
1480  lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
1481  lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
1482  lorybias_ = thePixelTemp_[index_id_].head.lorybias;
1483  lorxbias_ = thePixelTemp_[index_id_].head.lorxbias;
1484  if(flip_x) {lorxwidth_ = -lorxwidth_; lorxbias_ = -lorxbias_;}
1485  if(flip_y) {lorywidth_ = -lorywidth_; lorybias_ = -lorybias_;}
1486 
1487 
1488  }
1489 
1490  return success_;
1491 } // interpolate
1492 
1493 
1494 
1495 
1496 
1497 // ************************************************************************************************************
1503 // ************************************************************************************************************
1504 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta)
1505 {
1506  // Interpolate for a new set of track angles
1507 
1508  // Local variables
1509  float locBx = 1.f;
1510  if(cotbeta < 0.f) {locBx = -1.f;}
1511  float locBz = locBx;
1512  if(cotalpha < 0.f) {locBz = -locBx;}
1513  return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx);
1514 }
1515 
1516 
1517 
1518 // ************************************************************************************************************
1524 // ************************************************************************************************************
1525 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz)
1526 {
1527  // Interpolate for a new set of track angles
1528 
1529  // Local variables
1530  float locBx = 1.f;
1531  return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz, locBx);
1532 }
1533 
1534 
1535 // ************************************************************************************************************
1543 // ************************************************************************************************************
1544 void SiPixelTemplate::ysigma2(int fypix, int lypix, float sythr, float ysum[25], float ysig2[25])
1545 
1546 {
1547  // Interpolate using quantities already stored in the private variables
1548 
1549  // Local variables
1550  int i;
1551  float sigi, sigi2, sigi3, sigi4, symax, qscale, s25;
1552 
1553  // Make sure that input is OK
1554 
1555 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1556  if(fypix < 2 || fypix >= BYM2) {
1557  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with fypix = " << fypix << std::endl;
1558  }
1559 #else
1560  assert(fypix > 1 && fypix < BYM2);
1561 #endif
1562 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1563  if(lypix < fypix || lypix >= BYM2) {
1564  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with lypix/fypix = " << lypix << "/" << fypix << std::endl;
1565  }
1566 #else
1567  assert(lypix >= fypix && lypix < BYM2);
1568 #endif
1569 
1570  // Define the maximum signal to use in the parameterization
1571 
1572  symax = symax_;
1573  s25 = 0.5f*s50_;
1574  if(symax_ > syparmax_) {symax = syparmax_;}
1575 
1576  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1577 
1578  for(i=fypix-2; i<=lypix+2; ++i) {
1579  if(i < fypix || i > lypix) {
1580 
1581  // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
1582 
1583  ysig2[i] = s50_*s50_;
1584  } else {
1585  if(ysum[i] < symax) {
1586  sigi = ysum[i];
1587  qscale = 1.f;
1588  if(sigi < s25) sigi = s25;
1589  } else {
1590  sigi = symax;
1591  qscale = ysum[i]/symax;
1592  }
1593  sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
1594  if(i <= BHY) {
1595  ysig2[i] = (1.f-yratio_)*
1596  (yparl_[0][0]+yparl_[0][1]*sigi+yparl_[0][2]*sigi2+yparl_[0][3]*sigi3+yparl_[0][4]*sigi4)
1597  + yratio_*
1598  (yparh_[0][0]+yparh_[0][1]*sigi+yparh_[0][2]*sigi2+yparh_[0][3]*sigi3+yparh_[0][4]*sigi4);
1599  } else {
1600  ysig2[i] = (1.f-yratio_)*
1601  (yparl_[1][0]+yparl_[1][1]*sigi+yparl_[1][2]*sigi2+yparl_[1][3]*sigi3+yparl_[1][4]*sigi4)
1602  + yratio_*
1603  (yparh_[1][0]+yparh_[1][1]*sigi+yparh_[1][2]*sigi2+yparh_[1][3]*sigi3+yparh_[1][4]*sigi4);
1604  }
1605  ysig2[i] *=qscale;
1606  if(ysum[i] > sythr) {ysig2[i] = 1.e8;}
1607  if(ysig2[i] <= 0.f) {LOGERROR("SiPixelTemplate") << "neg y-error-squared, id = " << id_current_ << ", index = " << index_id_ <<
1608  ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_ << ", sigi = " << sigi << ENDL;}
1609  }
1610  }
1611 
1612  return;
1613 
1614 } // End ysigma2
1615 
1616 
1617 // ************************************************************************************************************
1623 // ************************************************************************************************************
1624 void SiPixelTemplate::ysigma2(float qpixel, int index, float& ysig2)
1625 
1626 {
1627  // Interpolate using quantities already stored in the private variables
1628 
1629  // Local variables
1630  float sigi, sigi2, sigi3, sigi4, symax, qscale, err2, s25;
1631 
1632  // Make sure that input is OK
1633 
1634 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1635  if(index < 2 || index >= BYM2) {
1636  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with index = " << index << std::endl;
1637  }
1638 #else
1639  assert(index > 1 && index < BYM2);
1640 #endif
1641 
1642  // Define the maximum signal to use in the parameterization
1643 
1644  symax = symax_;
1645  s25 = 0.5f*s50_;
1646  if(symax_ > syparmax_) {symax = syparmax_;}
1647 
1648  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1649 
1650  if(qpixel < symax) {
1651  sigi = qpixel;
1652  qscale = 1.f;
1653  if(sigi < s25) sigi = s25;
1654  } else {
1655  sigi = symax;
1656  qscale = qpixel/symax;
1657  }
1658  sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
1659  if(index <= BHY) {
1660  err2 = (1.f-yratio_)*
1661  (yparl_[0][0]+yparl_[0][1]*sigi+yparl_[0][2]*sigi2+yparl_[0][3]*sigi3+yparl_[0][4]*sigi4)
1662  + yratio_*
1663  (yparh_[0][0]+yparh_[0][1]*sigi+yparh_[0][2]*sigi2+yparh_[0][3]*sigi3+yparh_[0][4]*sigi4);
1664  } else {
1665  err2 = (1.f-yratio_)*
1666  (yparl_[1][0]+yparl_[1][1]*sigi+yparl_[1][2]*sigi2+yparl_[1][3]*sigi3+yparl_[1][4]*sigi4)
1667  + yratio_*
1668  (yparh_[1][0]+yparh_[1][1]*sigi+yparh_[1][2]*sigi2+yparh_[1][3]*sigi3+yparh_[1][4]*sigi4);
1669  }
1670  ysig2 =qscale*err2;
1671  if(ysig2 <= 0.f) {LOGERROR("SiPixelTemplate") << "neg y-error-squared, id = " << id_current_ << ", index = " << index_id_ <<
1672  ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_ << ", sigi = " << sigi << ENDL;}
1673 
1674  return;
1675 
1676 } // End ysigma2
1677 
1678 
1679 
1680 
1681 
1682 // ************************************************************************************************************
1690 // ************************************************************************************************************
1691 void SiPixelTemplate::xsigma2(int fxpix, int lxpix, float sxthr, float xsum[11], float xsig2[11])
1692 
1693 {
1694  // Interpolate using quantities already stored in the private variables
1695 
1696  // Local variables
1697  int i;
1698  float sigi, sigi2, sigi3, sigi4, yint, sxmax, x0, qscale, s25;
1699 
1700  // Make sure that input is OK
1701 
1702 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1703  if(fxpix < 2 || fxpix >= BXM2) {
1704  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xsigma2 called with fxpix = " << fxpix << std::endl;
1705  }
1706 #else
1707  assert(fxpix > 1 && fxpix < BXM2);
1708 #endif
1709 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1710  if(lxpix < fxpix || lxpix >= BXM2) {
1711  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xsigma2 called with lxpix/fxpix = " << lxpix << "/" << fxpix << std::endl;
1712  }
1713 #else
1714  assert(lxpix >= fxpix && lxpix < BXM2);
1715 #endif
1716 
1717  // Define the maximum signal to use in the parameterization
1718 
1719  sxmax = sxmax_;
1720  s25 = 0.5f*s50_;
1721  if(sxmax_ > sxparmax_) {sxmax = sxparmax_;}
1722 
1723  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1724 
1725  for(i=fxpix-2; i<=lxpix+2; ++i) {
1726  if(i < fxpix || i > lxpix) {
1727 
1728  // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
1729 
1730  xsig2[i] = s50_*s50_;
1731  } else {
1732  if(xsum[i] < sxmax) {
1733  sigi = xsum[i];
1734  qscale = 1.f;
1735  if(sigi < s25) sigi = s25;
1736  } else {
1737  sigi = sxmax;
1738  qscale = xsum[i]/sxmax;
1739  }
1740  sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
1741 
1742  // First, do the cotbeta interpolation
1743 
1744  if(i <= BHX) {
1745  yint = (1.f-yratio_)*
1746  (xparly0_[0][0]+xparly0_[0][1]*sigi+xparly0_[0][2]*sigi2+xparly0_[0][3]*sigi3+xparly0_[0][4]*sigi4)
1747  + yratio_*
1748  (xparhy0_[0][0]+xparhy0_[0][1]*sigi+xparhy0_[0][2]*sigi2+xparhy0_[0][3]*sigi3+xparhy0_[0][4]*sigi4);
1749  } else {
1750  yint = (1.f-yratio_)*
1751  (xparly0_[1][0]+xparly0_[1][1]*sigi+xparly0_[1][2]*sigi2+xparly0_[1][3]*sigi3+xparly0_[1][4]*sigi4)
1752  + yratio_*
1753  (xparhy0_[1][0]+xparhy0_[1][1]*sigi+xparhy0_[1][2]*sigi2+xparhy0_[1][3]*sigi3+xparhy0_[1][4]*sigi4);
1754  }
1755 
1756  // Next, do the cotalpha interpolation
1757 
1758  if(i <= BHX) {
1759  xsig2[i] = (1.f-xxratio_)*
1760  (xparl_[0][0]+xparl_[0][1]*sigi+xparl_[0][2]*sigi2+xparl_[0][3]*sigi3+xparl_[0][4]*sigi4)
1761  + xxratio_*
1762  (xparh_[0][0]+xparh_[0][1]*sigi+xparh_[0][2]*sigi2+xparh_[0][3]*sigi3+xparh_[0][4]*sigi4);
1763  } else {
1764  xsig2[i] = (1.f-xxratio_)*
1765  (xparl_[1][0]+xparl_[1][1]*sigi+xparl_[1][2]*sigi2+xparl_[1][3]*sigi3+xparl_[1][4]*sigi4)
1766  + xxratio_*
1767  (xparh_[1][0]+xparh_[1][1]*sigi+xparh_[1][2]*sigi2+xparh_[1][3]*sigi3+xparh_[1][4]*sigi4);
1768  }
1769 
1770  // Finally, get the mid-point value of the cotalpha function
1771 
1772  if(i <= BHX) {
1773  x0 = xpar0_[0][0]+xpar0_[0][1]*sigi+xpar0_[0][2]*sigi2+xpar0_[0][3]*sigi3+xpar0_[0][4]*sigi4;
1774  } else {
1775  x0 = xpar0_[1][0]+xpar0_[1][1]*sigi+xpar0_[1][2]*sigi2+xpar0_[1][3]*sigi3+xpar0_[1][4]*sigi4;
1776  }
1777 
1778  // Finally, rescale the yint value for cotalpha variation
1779 
1780  if(x0 != 0.f) {xsig2[i] = xsig2[i]/x0 * yint;}
1781  xsig2[i] *=qscale;
1782  if(xsum[i] > sxthr) {xsig2[i] = 1.e8;}
1783  if(xsig2[i] <= 0.f) {LOGERROR("SiPixelTemplate") << "neg x-error-squared, id = " << id_current_ << ", index = " << index_id_ <<
1784  ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_ << ", sigi = " << sigi << ENDL;}
1785  }
1786  }
1787 
1788  return;
1789 
1790 } // End xsigma2
1791 
1792 
1793 
1794 
1795 
1796 
1797 // ************************************************************************************************************
1801 // ************************************************************************************************************
1802 float SiPixelTemplate::yflcorr(int binq, float qfly)
1803 
1804 {
1805  // Interpolate using quantities already stored in the private variables
1806 
1807  // Local variables
1808  float qfl, qfl2, qfl3, qfl4, qfl5, dy;
1809 
1810  // Make sure that input is OK
1811 
1812 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1813  if(binq < 0 || binq > 3) {
1814  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yflcorr called with binq = " << binq << std::endl;
1815  }
1816 #else
1817  assert(binq >= 0 && binq < 4);
1818 #endif
1819 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1820  if(fabs((double)qfly) > 1.) {
1821  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yflcorr called with qfly = " << qfly << std::endl;
1822  }
1823 #else
1824  assert(fabs((double)qfly) <= 1.);
1825 #endif
1826 
1827  // Define the maximum signal to allow before de-weighting a pixel
1828 
1829  qfl = qfly;
1830 
1831  if(qfl < -0.9f) {qfl = -0.9f;}
1832  if(qfl > 0.9f) {qfl = 0.9f;}
1833 
1834  // Interpolate between the two polynomials
1835 
1836  qfl2 = qfl*qfl; qfl3 = qfl2*qfl; qfl4 = qfl3*qfl; qfl5 = qfl4*qfl;
1837  dy = (1.f-yratio_)*(yflparl_[binq][0]+yflparl_[binq][1]*qfl+yflparl_[binq][2]*qfl2+yflparl_[binq][3]*qfl3+yflparl_[binq][4]*qfl4+yflparl_[binq][5]*qfl5)
1838  + yratio_*(yflparh_[binq][0]+yflparh_[binq][1]*qfl+yflparh_[binq][2]*qfl2+yflparh_[binq][3]*qfl3+yflparh_[binq][4]*qfl4+yflparh_[binq][5]*qfl5);
1839 
1840  return dy;
1841 
1842 } // End yflcorr
1843 
1844 
1845 
1846 
1847 
1848 
1849 // ************************************************************************************************************
1853 // ************************************************************************************************************
1854 float SiPixelTemplate::xflcorr(int binq, float qflx)
1855 
1856 {
1857  // Interpolate using quantities already stored in the private variables
1858 
1859  // Local variables
1860  float qfl, qfl2, qfl3, qfl4, qfl5, dx;
1861 
1862  // Make sure that input is OK
1863 
1864 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1865  if(binq < 0 || binq > 3) {
1866  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xflcorr called with binq = " << binq << std::endl;
1867  }
1868 #else
1869  assert(binq >= 0 && binq < 4);
1870 #endif
1871 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1872  if(fabs((double)qflx) > 1.) {
1873  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xflcorr called with qflx = " << qflx << std::endl;
1874  }
1875 #else
1876  assert(fabs((double)qflx) <= 1.);
1877 #endif
1878 
1879  // Define the maximum signal to allow before de-weighting a pixel
1880 
1881  qfl = qflx;
1882 
1883  if(qfl < -0.9f) {qfl = -0.9f;}
1884  if(qfl > 0.9f) {qfl = 0.9f;}
1885 
1886  // Interpolate between the two polynomials
1887 
1888  qfl2 = qfl*qfl; qfl3 = qfl2*qfl; qfl4 = qfl3*qfl; qfl5 = qfl4*qfl;
1889  dx = (1.f - yxratio_)*((1.f-xxratio_)*(xflparll_[binq][0]+xflparll_[binq][1]*qfl+xflparll_[binq][2]*qfl2+xflparll_[binq][3]*qfl3+xflparll_[binq][4]*qfl4+xflparll_[binq][5]*qfl5)
1890  + xxratio_*(xflparlh_[binq][0]+xflparlh_[binq][1]*qfl+xflparlh_[binq][2]*qfl2+xflparlh_[binq][3]*qfl3+xflparlh_[binq][4]*qfl4+xflparlh_[binq][5]*qfl5))
1891  + yxratio_*((1.f-xxratio_)*(xflparhl_[binq][0]+xflparhl_[binq][1]*qfl+xflparhl_[binq][2]*qfl2+xflparhl_[binq][3]*qfl3+xflparhl_[binq][4]*qfl4+xflparhl_[binq][5]*qfl5)
1892  + xxratio_*(xflparhh_[binq][0]+xflparhh_[binq][1]*qfl+xflparhh_[binq][2]*qfl2+xflparhh_[binq][3]*qfl3+xflparhh_[binq][4]*qfl4+xflparhh_[binq][5]*qfl5));
1893 
1894  return dx;
1895 
1896 } // End xflcorr
1897 
1898 
1899 
1900 // ************************************************************************************************************
1905 // ************************************************************************************************************
1906 void SiPixelTemplate::ytemp(int fybin, int lybin, float ytemplate[41][BYSIZE])
1907 
1908 {
1909  // Retrieve already interpolated quantities
1910 
1911  // Local variables
1912  int i, j;
1913 
1914  // Verify that input parameters are in valid range
1915 
1916 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1917  if(fybin < 0 || fybin > 40) {
1918  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp called with fybin = " << fybin << std::endl;
1919  }
1920 #else
1921  assert(fybin >= 0 && fybin < 41);
1922 #endif
1923 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1924  if(lybin < 0 || lybin > 40) {
1925  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp called with lybin = " << lybin << std::endl;
1926  }
1927 #else
1928  assert(lybin >= 0 && lybin < 41);
1929 #endif
1930 
1931  // Build the y-template, the central 25 bins are here in all cases
1932 
1933  for(i=0; i<9; ++i) {
1934  for(j=0; j<BYSIZE; ++j) {
1935  ytemplate[i+16][j]=ytemp_[i][j];
1936  }
1937  }
1938  for(i=0; i<8; ++i) {
1939  ytemplate[i+8][BYM1] = 0.f;
1940  for(j=0; j<BYM1; ++j) {
1941  ytemplate[i+8][j]=ytemp_[i][j+1];
1942  }
1943  }
1944  for(i=1; i<9; ++i) {
1945  ytemplate[i+24][0] = 0.f;
1946  for(j=0; j<BYM1; ++j) {
1947  ytemplate[i+24][j+1]=ytemp_[i][j];
1948  }
1949  }
1950 
1951  // Add more bins if needed
1952 
1953  if(fybin < 8) {
1954  for(i=0; i<8; ++i) {
1955  ytemplate[i][BYM2] = 0.f;
1956  ytemplate[i][BYM1] = 0.f;
1957  for(j=0; j<BYM2; ++j) {
1958  ytemplate[i][j]=ytemp_[i][j+2];
1959  }
1960  }
1961  }
1962  if(lybin > 32) {
1963  for(i=1; i<9; ++i) {
1964  ytemplate[i+32][0] = 0.f;
1965  ytemplate[i+32][1] = 0.f;
1966  for(j=0; j<BYM2; ++j) {
1967  ytemplate[i+32][j+2]=ytemp_[i][j];
1968  }
1969  }
1970  }
1971 
1972  return;
1973 
1974 } // End ytemp
1975 
1976 
1977 
1978 // ************************************************************************************************************
1983 // ************************************************************************************************************
1984 void SiPixelTemplate::xtemp(int fxbin, int lxbin, float xtemplate[41][BXSIZE])
1985 
1986 {
1987  // Retrieve already interpolated quantities
1988 
1989  // Local variables
1990  int i, j;
1991 
1992  // Verify that input parameters are in valid range
1993 
1994 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1995  if(fxbin < 0 || fxbin > 40) {
1996  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp called with fxbin = " << fxbin << std::endl;
1997  }
1998 #else
1999  assert(fxbin >= 0 && fxbin < 41);
2000 #endif
2001 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2002  if(lxbin < 0 || lxbin > 40) {
2003  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp called with lxbin = " << lxbin << std::endl;
2004  }
2005 #else
2006  assert(lxbin >= 0 && lxbin < 41);
2007 #endif
2008 
2009  // Build the x-template, the central 25 bins are here in all cases
2010 
2011  for(i=0; i<9; ++i) {
2012  for(j=0; j<BXSIZE; ++j) {
2013  xtemplate[i+16][j]=xtemp_[i][j];
2014  }
2015  }
2016  for(i=0; i<8; ++i) {
2017  xtemplate[i+8][BXM1] = 0.f;
2018  for(j=0; j<BXM1; ++j) {
2019  xtemplate[i+8][j]=xtemp_[i][j+1];
2020  }
2021  }
2022  for(i=1; i<9; ++i) {
2023  xtemplate[i+24][0] = 0.f;
2024  for(j=0; j<BXM1; ++j) {
2025  xtemplate[i+24][j+1]=xtemp_[i][j];
2026  }
2027  }
2028 
2029  // Add more bins if needed
2030 
2031  if(fxbin < 8) {
2032  for(i=0; i<8; ++i) {
2033  xtemplate[i][BXM2] = 0.f;
2034  xtemplate[i][BXM1] = 0.f;
2035  for(j=0; j<BXM2; ++j) {
2036  xtemplate[i][j]=xtemp_[i][j+2];
2037  }
2038  }
2039  }
2040  if(lxbin > 32) {
2041  for(i=1; i<9; ++i) {
2042  xtemplate[i+32][0] = 0.f;
2043  xtemplate[i+32][1] = 0.f;
2044  for(j=0; j<BXM2; ++j) {
2045  xtemplate[i+32][j+2]=xtemp_[i][j];
2046  }
2047  }
2048  }
2049 
2050  return;
2051 
2052 } // End xtemp
2053 
2054 
2055 // ************************************************************************************************************
2057 // ************************************************************************************************************
2058 int SiPixelTemplate::cytemp()
2059 
2060 {
2061  // Retrieve already interpolated quantities
2062 
2063  // Local variables
2064  int j;
2065 
2066  // Analyze only pixels along the central entry
2067  // First, find the maximum signal and then work out to the edges
2068 
2069  float sigmax = 0.f;
2070  float qedge = 2.*s50_;
2071  int jmax = -1;
2072 
2073  for(j=0; j<BYSIZE; ++j) {
2074  if(ytemp_[4][j] > sigmax) {
2075  sigmax = ytemp_[4][j];
2076  jmax = j;
2077  }
2078  }
2079  if(sigmax < qedge) {qedge = s50_;}
2080  if(sigmax < qedge || jmax<1 || jmax>BYM2) {return -1;}
2081 
2082  // Now search forward and backward
2083 
2084  int jend = jmax;
2085 
2086  for(j=jmax+1; j<BYM1; ++j) {
2087  if(ytemp_[4][j] < qedge) break;
2088  jend = j;
2089  }
2090 
2091  int jbeg = jmax;
2092 
2093  for(j=jmax-1; j>0; --j) {
2094  if(ytemp_[4][j] < qedge) break;
2095  jbeg = j;
2096  }
2097 
2098  return (jbeg+jend)/2;
2099 
2100 } // End cytemp
2101 
2102 
2103 
2104 // ************************************************************************************************************
2106 // ************************************************************************************************************
2107 int SiPixelTemplate::cxtemp()
2108 
2109 {
2110  // Retrieve already interpolated quantities
2111 
2112  // Local variables
2113  int j;
2114 
2115  // Analyze only pixels along the central entry
2116  // First, find the maximum signal and then work out to the edges
2117 
2118  float sigmax = 0.f;
2119  float qedge = 2.*s50_;
2120  int jmax = -1;
2121 
2122  for(j=0; j<BXSIZE; ++j) {
2123  if(xtemp_[4][j] > sigmax) {
2124  sigmax = xtemp_[4][j];
2125  jmax = j;
2126  }
2127  }
2128  if(sigmax < qedge) {qedge = s50_;}
2129  if(sigmax < qedge || jmax<1 || jmax>BXM2) {return -1;}
2130 
2131  // Now search forward and backward
2132 
2133  int jend = jmax;
2134 
2135  for(j=jmax+1; j<BXM1; ++j) {
2136  if(xtemp_[4][j] < qedge) break;
2137  jend = j;
2138  }
2139 
2140  int jbeg = jmax;
2141 
2142  for(j=jmax-1; j>0; --j) {
2143  if(xtemp_[4][j] < qedge) break;
2144  jbeg = j;
2145  }
2146 
2147  return (jbeg+jend)/2;
2148 
2149 } // End cxtemp
2150 
2151 
2152 // ************************************************************************************************************
2156 // ************************************************************************************************************
2157 void SiPixelTemplate::ytemp3d_int(int nypix, int& nybins)
2158 
2159 {
2160 
2161  // Retrieve already interpolated quantities
2162 
2163  // Local variables
2164  int i, j, k;
2165  int ioff0, ioffp, ioffm;
2166 
2167  // Verify that input parameters are in valid range
2168 
2169 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2170  if(nypix < 1 || nypix >= BYM3) {
2171  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp3d called with nypix = " << nypix << std::endl;
2172  }
2173 #else
2174  assert(nypix > 0 && nypix < BYM3);
2175 #endif
2176 
2177  // Calculate the size of the shift in pixels needed to span the entire cluster
2178 
2179  float diff = fabsf(nypix - clsleny_)/2. + 1.f;
2180  int nshift = (int)diff;
2181  if((diff - nshift) > 0.5f) {++nshift;}
2182 
2183  // Calculate the number of bins needed to specify each hit range
2184 
2185  nybins_ = 9 + 16*nshift;
2186 
2187  // Create a 2-d working template with the correct size
2188 
2189  temp2dy_.resize(boost::extents[nybins_][BYSIZE]);
2190 
2191  // The 9 central bins are copied from the interpolated private store
2192 
2193  ioff0 = 8*nshift;
2194 
2195  for(i=0; i<9; ++i) {
2196  for(j=0; j<BYSIZE; ++j) {
2197  temp2dy_[i+ioff0][j]=ytemp_[i][j];
2198  }
2199  }
2200 
2201  // Add the +- shifted templates
2202 
2203  for(k=1; k<=nshift; ++k) {
2204  ioffm=ioff0-k*8;
2205  for(i=0; i<8; ++i) {
2206  for(j=0; j<k; ++j) {
2207  temp2dy_[i+ioffm][BYM1-j] = 0.f;
2208  }
2209  for(j=0; j<BYSIZE-k; ++j) {
2210  temp2dy_[i+ioffm][j]=ytemp_[i][j+k];
2211  }
2212  }
2213  ioffp=ioff0+k*8;
2214  for(i=1; i<9; ++i) {
2215  for(j=0; j<k; ++j) {
2216  temp2dy_[i+ioffp][j] = 0.f;
2217  }
2218  for(j=0; j<BYSIZE-k; ++j) {
2219  temp2dy_[i+ioffp][j+k]=ytemp_[i][j];
2220  }
2221  }
2222  }
2223 
2224  nybins = nybins_;
2225  return;
2226 
2227 } // End ytemp3d_int
2228 
2229 
2230 // ************************************************************************************************************
2234 // ************************************************************************************************************
2235 void SiPixelTemplate::ytemp3d(int i, int j, std::vector<float>& ytemplate)
2236 
2237 {
2238  // Sum two 2-d templates to make the 3-d template
2239  if(i >= 0 && i < nybins_ && j <= i) {
2240  for(int k=0; k<BYSIZE; ++k) {
2241  ytemplate[k]=temp2dy_[i][k]+temp2dy_[j][k];
2242  }
2243  } else {
2244  for(int k=0; k<BYSIZE; ++k) {
2245  ytemplate[k]=0.;
2246  }
2247  }
2248 
2249  return;
2250 
2251 } // End ytemp3d
2252 
2253 
2254 // ************************************************************************************************************
2258 // ************************************************************************************************************
2259 void SiPixelTemplate::xtemp3d_int(int nxpix, int& nxbins)
2260 
2261 {
2262  // Retrieve already interpolated quantities
2263 
2264  // Local variables
2265  int i, j, k;
2266  int ioff0, ioffp, ioffm;
2267 
2268  // Verify that input parameters are in valid range
2269 
2270 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2271  if(nxpix < 1 || nxpix >= BXM3) {
2272  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp3d called with nxpix = " << nxpix << std::endl;
2273  }
2274 #else
2275  assert(nxpix > 0 && nxpix < BXM3);
2276 #endif
2277 
2278  // Calculate the size of the shift in pixels needed to span the entire cluster
2279 
2280  float diff = fabsf(nxpix - clslenx_)/2.f + 1.f;
2281  int nshift = (int)diff;
2282  if((diff - nshift) > 0.5f) {++nshift;}
2283 
2284  // Calculate the number of bins needed to specify each hit range
2285 
2286  nxbins_ = 9 + 16*nshift;
2287 
2288  // Create a 2-d working template with the correct size
2289 
2290  temp2dx_.resize(boost::extents[nxbins_][BXSIZE]);
2291 
2292  // The 9 central bins are copied from the interpolated private store
2293 
2294  ioff0 = 8*nshift;
2295 
2296  for(i=0; i<9; ++i) {
2297  for(j=0; j<BXSIZE; ++j) {
2298  temp2dx_[i+ioff0][j]=xtemp_[i][j];
2299  }
2300  }
2301 
2302  // Add the +- shifted templates
2303 
2304  for(k=1; k<=nshift; ++k) {
2305  ioffm=ioff0-k*8;
2306  for(i=0; i<8; ++i) {
2307  for(j=0; j<k; ++j) {
2308  temp2dx_[i+ioffm][BXM1-j] = 0.f;
2309  }
2310  for(j=0; j<BXSIZE-k; ++j) {
2311  temp2dx_[i+ioffm][j]=xtemp_[i][j+k];
2312  }
2313  }
2314  ioffp=ioff0+k*8;
2315  for(i=1; i<9; ++i) {
2316  for(j=0; j<k; ++j) {
2317  temp2dx_[i+ioffp][j] = 0.f;
2318  }
2319  for(j=0; j<BXSIZE-k; ++j) {
2320  temp2dx_[i+ioffp][j+k]=xtemp_[i][j];
2321  }
2322  }
2323  }
2324 
2325  nxbins = nxbins_;
2326 
2327  return;
2328 
2329 } // End xtemp3d_int
2330 
2331 
2332 
2333 // ************************************************************************************************************
2337 // ************************************************************************************************************
2338 void SiPixelTemplate::xtemp3d(int i, int j, std::vector<float>& xtemplate)
2339 
2340 {
2341  // Sum two 2-d templates to make the 3-d template
2342  if(i >= 0 && i < nxbins_ && j <= i) {
2343  for(int k=0; k<BXSIZE; ++k) {
2344  xtemplate[k]=temp2dx_[i][k]+temp2dx_[j][k];
2345  }
2346  } else {
2347  for(int k=0; k<BXSIZE; ++k) {
2348  xtemplate[k]=0.;
2349  }
2350  }
2351 
2352  return;
2353 
2354 } // End xtemp3d
2355 
2356 
2357 // ************************************************************************************************************
2385 // ************************************************************************************************************
2386 int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float locBz, float locBx, float qclus, float& pixmx,
2387  float& sigmay, float& deltay, float& sigmax, float& deltax,float& sy1, float& dy1, float& sy2,
2388  float& dy2, float& sx1, float& dx1, float& sx2, float& dx2) // float& lorywidth, float& lorxwidth)
2389 {
2390  // Interpolate for a new set of track angles
2391 
2392 
2393  // Find the index corresponding to id
2394 
2395  int index = -1;
2396  for( int i=0; i<(int)thePixelTemp_.size(); ++i) {
2397  if(id == thePixelTemp_[i].head.ID) {
2398  index = i;
2399  break;
2400  }
2401  }
2402 
2403 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2404  if(index < 0 || index >= (int)thePixelTemp_.size()) {
2405  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::qbin can't find needed template ID = " << id << std::endl;
2406  }
2407 #else
2408  assert(index >= 0 && index < (int)thePixelTemp_.size());
2409 #endif
2410 
2411  //
2412 
2413 
2414  auto const & templ = thePixelTemp_[index];
2415 
2416  // Interpolate the absolute value of cot(beta)
2417 
2418  auto acotb = std::abs(cotbeta);
2419 
2420  // qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
2421 
2422  auto cotalpha0 = thePixelTemp_[index].enty[0].cotalpha;
2423  auto qcorrect=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta*cotbeta+cotalpha0*cotalpha0));
2424 
2425  // for some cosmics, the ususal gymnastics are incorrect
2426 
2427  float cota = cotalpha;
2428  float cotb = abs_cotb_;
2429  bool flip_x = false;
2430  bool flip_y = false;
2431  switch(thePixelTemp_[index_id_].head.Dtype) {
2432  case 0:
2433  if(cotbeta < 0.f) {flip_y = true;}
2434  break;
2435  case 1:
2436  if(locBz < 0.f) {
2437  cotb = cotbeta;
2438  } else {
2439  cotb = -cotbeta;
2440  flip_y = true;
2441  }
2442  break;
2443  case 2:
2444  case 3:
2445  case 4:
2446  case 5:
2447  if(locBx*locBz < 0.f) {
2448  cota = -cotalpha;
2449  flip_x = true;
2450  }
2451  if(locBx > 0.f) {
2452  cotb = cotbeta;
2453  } else {
2454  cotb = -cotbeta;
2455  flip_y = true;
2456  }
2457  break;
2458  default:
2459 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2460  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
2461 #else
2462  std::cout << "SiPixelTemplate::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
2463 #endif
2464  }
2465 
2466 
2467  // Copy the charge scaling factor to the private variable
2468 
2469 
2470  auto qscale = thePixelTemp_[index].head.qscale;
2471 
2472 
2473  /*
2474  lorywidth = thePixelTemp_[index].head.lorywidth;
2475  if(locBz > 0.f) {lorywidth = -lorywidth;}
2476  lorxwidth = thePixelTemp_[index].head.lorxwidth;
2477  */
2478 
2479 
2480  auto Ny = thePixelTemp_[index].head.NTy;
2481  auto Nyx = thePixelTemp_[index].head.NTyx;
2482  auto Nxx = thePixelTemp_[index].head.NTxx;
2483 
2484 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2485  if(Ny < 2 || Nyx < 1 || Nxx < 2) {
2486  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
2487  }
2488 #else
2489  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
2490 #endif
2491 
2492  // next, loop over all y-angle entries
2493 
2494  auto ilow = 0;
2495  auto ihigh = 0;
2496  auto yratio = 0.f;
2497 
2498  {
2499  auto j = std::lower_bound(templ.cotbetaY,templ.cotbetaY+Ny,cotb);
2500  if (j==templ.cotbetaY+Ny) { --j; yratio = 1.f; }
2501  else if (j==templ.cotbetaY) { ++j; yratio = 0.f;}
2502  else { yratio = (cotb - (*(j-1)) )/( (*j) - (*(j-1)) ) ; }
2503 
2504  ihigh = j-templ.cotbetaY;
2505  ilow = ihigh-1;
2506  }
2507 
2508 
2509 
2510  // Interpolate/store all y-related quantities (flip displacements when flip_y)
2511 
2512  dy1 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].dyone + yratio*thePixelTemp_[index].enty[ihigh].dyone;
2513  if(flip_y) {dy1 = -dy1;}
2514  sy1 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].syone + yratio*thePixelTemp_[index].enty[ihigh].syone;
2515  dy2 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].dytwo + yratio*thePixelTemp_[index].enty[ihigh].dytwo;
2516  if(flip_y) {dy2 = -dy2;}
2517  sy2 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].sytwo + yratio*thePixelTemp_[index].enty[ihigh].sytwo;
2518 
2519  auto qavg = (1.f - yratio)*thePixelTemp_[index].enty[ilow].qavg + yratio*thePixelTemp_[index].enty[ihigh].qavg;
2520  qavg *= qcorrect;
2521  auto qmin = (1.f - yratio)*thePixelTemp_[index].enty[ilow].qmin + yratio*thePixelTemp_[index].enty[ihigh].qmin;
2522  qmin *= qcorrect;
2523  auto qmin2 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].qmin2 + yratio*thePixelTemp_[index].enty[ihigh].qmin2;
2524  qmin2 *= qcorrect;
2525 
2526 
2527 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2528  if(qavg <= 0.f || qmin <= 0.f) {
2529  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::qbin, qavg or qmin <= 0,"
2530  << " Probably someone called the generic pixel reconstruction with an illegal trajectory state" << std::endl;
2531  }
2532 #else
2533  assert(qavg > 0.f && qmin > 0.f);
2534 #endif
2535 
2536  // Scale the input charge to account for differences between pixelav and CMSSW simulation or data
2537 
2538  auto qtotal = qscale*qclus;
2539 
2540  // uncertainty and final corrections depend upon total charge bin
2541  auto fq = qtotal/qavg;
2542  int binq;
2543 
2544  if(fq > fbin_[0]) {
2545  binq=0;
2546  //std::cout<<" fq "<<fq<<" "<<qtotal<<" "<<qavg<<" "<<qclus<<" "<<qscale<<" "
2547  // <<fbin_[0]<<" "<<fbin_[1]<<" "<<fbin_[2]<<std::endl;
2548  } else {
2549  if(fq > fbin_[1]) {
2550  binq=1;
2551  } else {
2552  if(fq > fbin_[2]) {
2553  binq=2;
2554  } else {
2555  binq=3;
2556  }
2557  }
2558  }
2559 
2560  auto yavggen =(1.f - yratio)*thePixelTemp_[index].enty[ilow].yavggen[binq] + yratio*thePixelTemp_[index].enty[ihigh].yavggen[binq];
2561  if(flip_y) {yavggen = -yavggen;}
2562  auto yrmsgen =(1.f - yratio)*thePixelTemp_[index].enty[ilow].yrmsgen[binq] + yratio*thePixelTemp_[index].enty[ihigh].yrmsgen[binq];
2563 
2564 
2565  // next, loop over all x-angle entries, first, find relevant y-slices
2566 
2567  auto iylow = 0;
2568  auto iyhigh = 0;
2569  auto yxratio = 0.f;
2570 
2571 
2572  {
2573  auto j = std::lower_bound(templ.cotbetaX,templ.cotbetaX+Nyx,acotb);
2574  if (j==templ.cotbetaX+Nyx) { --j; yxratio = 1.f; }
2575  else if (j==templ.cotbetaX) { ++j; yxratio = 0.f;}
2576  else { yxratio = (acotb - (*(j-1)) )/( (*j) - (*(j-1)) ) ; }
2577 
2578  iyhigh = j-templ.cotbetaX;
2579  iylow = iyhigh -1;
2580  }
2581 
2582 
2583 
2584  // Unneded: ilow = ihigh = 0;
2585  auto xxratio = 0.f;
2586 
2587  {
2588  auto j = std::lower_bound(templ.cotalphaX,templ.cotalphaX+Nxx,cota);
2589  if (j==templ.cotalphaX+Nxx) { --j; xxratio = 1.f; }
2590  else if (j==templ.cotalphaX) { ++j; xxratio = 0.f;}
2591  else { xxratio = (cota - (*(j-1)) )/( (*j) - (*(j-1)) ) ; }
2592 
2593  ihigh = j-templ.cotalphaX;
2594  ilow = ihigh-1;
2595  }
2596 
2597 
2598 
2599  dx1 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].dxone + xxratio*thePixelTemp_[index].entx[0][ihigh].dxone;
2600  if(flip_x) {dx1 = -dx1;}
2601  sx1 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].sxone + xxratio*thePixelTemp_[index].entx[0][ihigh].sxone;
2602  dx2 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].dxtwo + xxratio*thePixelTemp_[index].entx[0][ihigh].dxtwo;
2603  if(flip_x) {dx2 = -dx2;}
2604  sx2 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].sxtwo + xxratio*thePixelTemp_[index].entx[0][ihigh].sxtwo;
2605 
2606  // pixmax is the maximum allowed pixel charge (used for truncation)
2607 
2608  pixmx=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].pixmax + xxratio*thePixelTemp_[index].entx[iylow][ihigh].pixmax)
2609  +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].pixmax + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].pixmax);
2610 
2611  auto xavggen = (1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].xavggen[binq] + xxratio*thePixelTemp_[index].entx[iylow][ihigh].xavggen[binq])
2612  +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].xavggen[binq] + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].xavggen[binq]);
2613  if(flip_x) {xavggen = -xavggen;}
2614 
2615  auto xrmsgen = (1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].xrmsgen[binq] + xxratio*thePixelTemp_[index].entx[iylow][ihigh].xrmsgen[binq])
2616  +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].xrmsgen[binq] + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].xrmsgen[binq]);
2617 
2618 
2619 
2620  // Take the errors and bias from the correct charge bin
2621 
2622  sigmay = yrmsgen; deltay = yavggen;
2623 
2624  sigmax = xrmsgen; deltax = xavggen;
2625 
2626  // If the charge is too small (then flag it)
2627 
2628  if(qtotal < 0.95f*qmin) {binq = 5;} else {if(qtotal < 0.95f*qmin2) {binq = 4;}}
2629 
2630  return binq;
2631 
2632 } // qbin
2633 
2634 
2635 // ************************************************************************************************************
2663 // ************************************************************************************************************
2664 int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float& pixmx,
2665  float& sigmay, float& deltay, float& sigmax, float& deltax,float& sy1, float& dy1, float& sy2,
2666  float& dy2, float& sx1, float& dx1, float& sx2, float& dx2) // float& lorywidth, float& lorxwidth)
2667 {
2668  // Interpolate for a new set of track angles
2669 
2670  // Local variables
2671  float locBx = 1.f; // lorywidth, lorxwidth;
2672 
2673  return SiPixelTemplate::qbin(id, cotalpha, cotbeta, locBz, locBx, qclus, pixmx, sigmay, deltay, sigmax, deltax,
2674  sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2); // , lorywidth, lorxwidth);
2675 
2676 } // qbin
2677 
2678 // ************************************************************************************************************
2700 // ************************************************************************************************************
2701 /*
2702  int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float& pixmx, float& sigmay, float& deltay, float& sigmax, float& deltax,
2703  float& sy1, float& dy1, float& sy2, float& dy2, float& sx1, float& dx1, float& sx2, float& dx2)
2704 
2705  {
2706  float lorywidth, lorxwidth;
2707  return SiPixelTemplate::qbin(id, cotalpha, cotbeta, locBz, qclus, pixmx, sigmay, deltay, sigmax, deltax,
2708  sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, lorywidth, lorxwidth);
2709 
2710  } // qbin
2711  */
2712 
2713 // ************************************************************************************************************
2720 // ************************************************************************************************************
2721 int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float qclus)
2722 {
2723  // Interpolate for a new set of track angles
2724 
2725  // Local variables
2726  float pixmx, sigmay, deltay, sigmax, deltax, sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, locBz; // lorywidth, lorxwidth;
2727  // Local variables
2728  float locBx = 1.f;
2729  if(cotbeta < 0.f) {locBx = -1.f;}
2730  locBz = locBx;
2731  if(cotalpha < 0.f) {locBz = -locBx;}
2732 
2733  return SiPixelTemplate::qbin(id, cotalpha, cotbeta, locBz, locBx, qclus, pixmx, sigmay, deltay, sigmax, deltax,
2734  sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2); // , lorywidth, lorxwidth);
2735 
2736 } // qbin
2737 
2738 // ************************************************************************************************************
2744 // ************************************************************************************************************
2745 int SiPixelTemplate::qbin(int id, float cotbeta, float qclus)
2746 {
2747  // Interpolate for a new set of track angles
2748 
2749  // Local variables
2750  float pixmx, sigmay, deltay, sigmax, deltax, sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, locBz, locBx; //, lorywidth, lorxwidth;
2751  const float cotalpha = 0.f;
2752  locBx = 1.f;
2753  if(cotbeta < 0.f) {locBx = -1.f;}
2754  locBz = locBx;
2755  if(cotalpha < 0.f) {locBz = -locBx;}
2756  return SiPixelTemplate::qbin(id, cotalpha, cotbeta, locBz, locBx, qclus, pixmx, sigmay, deltay, sigmax, deltax,
2757  sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2); // , lorywidth, lorxwidth);
2758 
2759 } // qbin
2760 
2761 
2762 
2763 
2764 // ************************************************************************************************************
2776 // ************************************************************************************************************
2777 void SiPixelTemplate::temperrors(int id, float cotalpha, float cotbeta, int qBin, float& sigmay, float& sigmax, float& sy1, float& sy2, float& sx1, float& sx2)
2778 
2779 {
2780  // Interpolate for a new set of track angles
2781 
2782  // Local variables
2783  int i;
2784  int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, index;
2785  float yratio, yxratio, xxratio;
2786  float acotb, cotb;
2787  float yrms, xrms;
2788  //bool flip_y;
2789 
2790 
2791  // Find the index corresponding to id
2792 
2793  index = -1;
2794  for(i=0; i<(int)thePixelTemp_.size(); ++i) {
2795 
2796  if(id == thePixelTemp_[i].head.ID) {
2797 
2798  index = i;
2799  break;
2800  }
2801  }
2802 
2803 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2804  if(index < 0 || index >= (int)thePixelTemp_.size()) {
2805  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors can't find needed template ID = " << id << std::endl;
2806  }
2807 #else
2808  assert(index >= 0 && index < (int)thePixelTemp_.size());
2809 #endif
2810 
2811 
2812 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2813  if(qBin < 0 || qBin > 5) {
2814  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors called with illegal qBin = " << qBin << std::endl;
2815  }
2816 #else
2817  assert(qBin >= 0 && qBin < 6);
2818 #endif
2819 
2820  // The error information for qBin > 3 is taken to be the same as qBin=3
2821 
2822  if(qBin > 3) {qBin = 3;}
2823  //
2824 
2825  // Interpolate the absolute value of cot(beta)
2826 
2827  acotb = fabs((double)cotbeta);
2828  // cotb = cotbeta; // &&& check with Morris, we reassign it below.
2829 
2830  // for some cosmics, the ususal gymnastics are incorrect
2831 
2832  // if(thePixelTemp_[index].head.Dtype == 0) {
2833  cotb = acotb;
2834  //flip_y = false;
2835  //if(cotbeta < 0.f) {flip_y = true;}
2836  // } else {
2837  // if(locBz < 0.f) {
2838  // cotb = cotbeta;
2839  // flip_y = false;
2840  // } else {
2841  // cotb = -cotbeta;
2842  // flip_y = true;
2843  // }
2844  // }
2845 
2846  // Copy the charge scaling factor to the private variable
2847 
2848  Ny = thePixelTemp_[index].head.NTy;
2849  Nyx = thePixelTemp_[index].head.NTyx;
2850  Nxx = thePixelTemp_[index].head.NTxx;
2851 
2852 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2853  if(Ny < 2 || Nyx < 1 || Nxx < 2) {
2854  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
2855  }
2856 #else
2857  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
2858 #endif
2859 
2860  // next, loop over all y-angle entries
2861 
2862  ilow = 0;
2863  yratio = 0.f;
2864 
2865  if(cotb >= thePixelTemp_[index].enty[Ny-1].cotbeta) {
2866 
2867  ilow = Ny-2;
2868  yratio = 1.f;
2869 
2870  } else {
2871 
2872  if(cotb >= thePixelTemp_[index].enty[0].cotbeta) {
2873 
2874  for (i=0; i<Ny-1; ++i) {
2875 
2876  if( thePixelTemp_[index].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index].enty[i+1].cotbeta) {
2877 
2878  ilow = i;
2879  yratio = (cotb - thePixelTemp_[index].enty[i].cotbeta)/(thePixelTemp_[index].enty[i+1].cotbeta - thePixelTemp_[index].enty[i].cotbeta);
2880  break;
2881  }
2882  }
2883  }
2884  }
2885 
2886  ihigh=ilow + 1;
2887 
2888  // Interpolate/store all y-related quantities (flip displacements when flip_y)
2889 
2890  sy1 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].syone + yratio*thePixelTemp_[index].enty[ihigh].syone;
2891  sy2 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].sytwo + yratio*thePixelTemp_[index].enty[ihigh].sytwo;
2892  yrms=(1.f - yratio)*thePixelTemp_[index].enty[ilow].yrms[qBin] + yratio*thePixelTemp_[index].enty[ihigh].yrms[qBin];
2893 
2894 
2895  // next, loop over all x-angle entries, first, find relevant y-slices
2896 
2897  iylow = 0;
2898  yxratio = 0.f;
2899 
2900  if(acotb >= thePixelTemp_[index].entx[Nyx-1][0].cotbeta) {
2901 
2902  iylow = Nyx-2;
2903  yxratio = 1.f;
2904 
2905  } else if(acotb >= thePixelTemp_[index].entx[0][0].cotbeta) {
2906 
2907  for (i=0; i<Nyx-1; ++i) {
2908 
2909  if( thePixelTemp_[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp_[index].entx[i+1][0].cotbeta) {
2910 
2911  iylow = i;
2912  yxratio = (acotb - thePixelTemp_[index].entx[i][0].cotbeta)/(thePixelTemp_[index].entx[i+1][0].cotbeta - thePixelTemp_[index].entx[i][0].cotbeta);
2913  break;
2914  }
2915  }
2916  }
2917 
2918  iyhigh=iylow + 1;
2919 
2920  ilow = 0;
2921  xxratio = 0.f;
2922 
2923  if(cotalpha >= thePixelTemp_[index].entx[0][Nxx-1].cotalpha) {
2924 
2925  ilow = Nxx-2;
2926  xxratio = 1.f;
2927 
2928  } else {
2929 
2930  if(cotalpha >= thePixelTemp_[index].entx[0][0].cotalpha) {
2931 
2932  for (i=0; i<Nxx-1; ++i) {
2933 
2934  if( thePixelTemp_[index].entx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp_[index].entx[0][i+1].cotalpha) {
2935 
2936  ilow = i;
2937  xxratio = (cotalpha - thePixelTemp_[index].entx[0][i].cotalpha)/(thePixelTemp_[index].entx[0][i+1].cotalpha - thePixelTemp_[index].entx[0][i].cotalpha);
2938  break;
2939  }
2940  }
2941  }
2942  }
2943 
2944  ihigh=ilow + 1;
2945 
2946  sx1 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].sxone + xxratio*thePixelTemp_[index].entx[0][ihigh].sxone;
2947  sx2 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].sxtwo + xxratio*thePixelTemp_[index].entx[0][ihigh].sxtwo;
2948 
2949  xrms=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].xrms[qBin] + xxratio*thePixelTemp_[index].entx[iylow][ihigh].xrms[qBin])
2950  +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].xrms[qBin] + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].xrms[qBin]);
2951 
2952 
2953 
2954 
2955  // Take the errors and bias from the correct charge bin
2956 
2957  sigmay = yrms;
2958 
2959  sigmax = xrms;
2960 
2961  return;
2962 
2963 } // temperrors
2964 
2965 // ************************************************************************************************************
2975 // ************************************************************************************************************
2976 void SiPixelTemplate::qbin_dist(int id, float cotalpha, float cotbeta, float qbin_frac[4], float& ny1_frac, float& ny2_frac, float& nx1_frac, float& nx2_frac)
2977 
2978 {
2979  // Interpolate for a new set of track angles
2980 
2981  // Local variables
2982  int i;
2983  int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, index;
2984  float yratio, yxratio, xxratio;
2985  float acotb, cotb;
2986  float qfrac[4];
2987  //bool flip_y;
2988 
2989  // Find the index corresponding to id
2990 
2991  index = -1;
2992  for(i=0; i<(int)thePixelTemp_.size(); ++i) {
2993 
2994  if(id == thePixelTemp_[i].head.ID) {
2995 
2996  index = i;
2997  // id_current_ = id;
2998  break;
2999  }
3000  }
3001 
3002 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3003  if(index < 0 || index >= (int)thePixelTemp_.size()) {
3004  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors can't find needed template ID = " << id << std::endl;
3005  }
3006 #else
3007  assert(index >= 0 && index < (int)thePixelTemp_.size());
3008 #endif
3009 
3010  //
3011 
3012  // Interpolate the absolute value of cot(beta)
3013 
3014  acotb = fabs((double)cotbeta);
3015  // cotb = cotbeta; // &&& check with Morris, we reassign it below.
3016 
3017 
3018  // for some cosmics, the ususal gymnastics are incorrect
3019 
3020  // if(thePixelTemp_[index].head.Dtype == 0) {
3021  cotb = acotb;
3022  //flip_y = false;
3023  //if(cotbeta < 0.f) {flip_y = true;}
3024  // } else {
3025  // if(locBz < 0.f) {
3026  // cotb = cotbeta;
3027  // flip_y = false;
3028  // } else {
3029  // cotb = -cotbeta;
3030  // flip_y = true;
3031  // }
3032  // }
3033 
3034  // Copy the charge scaling factor to the private variable
3035 
3036  Ny = thePixelTemp_[index].head.NTy;
3037  Nyx = thePixelTemp_[index].head.NTyx;
3038  Nxx = thePixelTemp_[index].head.NTxx;
3039 
3040 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3041  if(Ny < 2 || Nyx < 1 || Nxx < 2) {
3042  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
3043  }
3044 #else
3045  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
3046 #endif
3047 
3048  // next, loop over all y-angle entries
3049 
3050  ilow = 0;
3051  yratio = 0.f;
3052 
3053  if(cotb >= thePixelTemp_[index].enty[Ny-1].cotbeta) {
3054 
3055  ilow = Ny-2;
3056  yratio = 1.f;
3057 
3058  } else {
3059 
3060  if(cotb >= thePixelTemp_[index].enty[0].cotbeta) {
3061 
3062  for (i=0; i<Ny-1; ++i) {
3063 
3064  if( thePixelTemp_[index].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index].enty[i+1].cotbeta) {
3065 
3066  ilow = i;
3067  yratio = (cotb - thePixelTemp_[index].enty[i].cotbeta)/(thePixelTemp_[index].enty[i+1].cotbeta - thePixelTemp_[index].enty[i].cotbeta);
3068  break;
3069  }
3070  }
3071  }
3072  }
3073 
3074  ihigh=ilow + 1;
3075 
3076  // Interpolate/store all y-related quantities (flip displacements when flip_y)
3077  ny1_frac = (1.f - yratio)*thePixelTemp_[index].enty[ilow].fracyone + yratio*thePixelTemp_[index].enty[ihigh].fracyone;
3078  ny2_frac = (1.f - yratio)*thePixelTemp_[index].enty[ilow].fracytwo + yratio*thePixelTemp_[index].enty[ihigh].fracytwo;
3079 
3080  // next, loop over all x-angle entries, first, find relevant y-slices
3081 
3082  iylow = 0;
3083  yxratio = 0.f;
3084 
3085  if(acotb >= thePixelTemp_[index].entx[Nyx-1][0].cotbeta) {
3086 
3087  iylow = Nyx-2;
3088  yxratio = 1.f;
3089 
3090  } else if(acotb >= thePixelTemp_[index].entx[0][0].cotbeta) {
3091 
3092  for (i=0; i<Nyx-1; ++i) {
3093 
3094  if( thePixelTemp_[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp_[index].entx[i+1][0].cotbeta) {
3095 
3096  iylow = i;
3097  yxratio = (acotb - thePixelTemp_[index].entx[i][0].cotbeta)/(thePixelTemp_[index].entx[i+1][0].cotbeta - thePixelTemp_[index].entx[i][0].cotbeta);
3098  break;
3099  }
3100  }
3101  }
3102 
3103  iyhigh=iylow + 1;
3104 
3105  ilow = 0;
3106  xxratio = 0.f;
3107 
3108  if(cotalpha >= thePixelTemp_[index].entx[0][Nxx-1].cotalpha) {
3109 
3110  ilow = Nxx-2;
3111  xxratio = 1.f;
3112 
3113  } else {
3114 
3115  if(cotalpha >= thePixelTemp_[index].entx[0][0].cotalpha) {
3116 
3117  for (i=0; i<Nxx-1; ++i) {
3118 
3119  if( thePixelTemp_[index].entx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp_[index].entx[0][i+1].cotalpha) {
3120 
3121  ilow = i;
3122  xxratio = (cotalpha - thePixelTemp_[index].entx[0][i].cotalpha)/(thePixelTemp_[index].entx[0][i+1].cotalpha - thePixelTemp_[index].entx[0][i].cotalpha);
3123  break;
3124  }
3125  }
3126  }
3127  }
3128 
3129  ihigh=ilow + 1;
3130 
3131  for(i=0; i<3; ++i) {
3132  qfrac[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].qbfrac[i] + xxratio*thePixelTemp_[index].entx[iylow][ihigh].qbfrac[i])
3133  +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].qbfrac[i] + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].qbfrac[i]);
3134  }
3135  nx1_frac = (1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].fracxone + xxratio*thePixelTemp_[index].entx[iylow][ihigh].fracxone)
3136  +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].fracxone + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].fracxone);
3137  nx2_frac = (1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].fracxtwo + xxratio*thePixelTemp_[index].entx[iylow][ihigh].fracxtwo)
3138  +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].fracxtwo + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].fracxtwo);
3139 
3140 
3141 
3142  qbin_frac[0] = qfrac[0];
3143  qbin_frac[1] = qbin_frac[0] + qfrac[1];
3144  qbin_frac[2] = qbin_frac[1] + qfrac[2];
3145  qbin_frac[3] = 1.f;
3146  return;
3147 
3148 } // qbin
3149 
3150 // *************************************************************************************************************************************
3152 
3158 // *************************************************************************************************************************************
3159 
3160 bool SiPixelTemplate::simpletemplate2D(float xhit, float yhit, std::vector<bool>& ydouble, std::vector<bool>& xdouble, float template2d[BXM2][BYM2])
3161 {
3162 
3163  // Local variables
3164 
3165  float x0, y0, xf, yf, xi, yi, sf, si, s0, qpix, slopey, slopex, ds;
3166  int i, j, jpix0, ipix0, jpixf, ipixf, jpix, ipix, nx, ny, anx, any, jmax, imax;
3167  float qtotal;
3168  // double path;
3169  std::list<SimplePixel> list;
3170  std::list<SimplePixel>::iterator listIter, listEnd;
3171 
3172  // Calculate the entry and exit points for the line charge from the track
3173 
3174  x0 = xhit - 0.5*zsize_*cota_current_;
3175  y0 = yhit - 0.5*zsize_*cotb_current_;
3176 
3177  jpix0 = floor(x0/xsize_)+1;
3178  ipix0 = floor(y0/ysize_)+1;
3179 
3180  if(jpix0 < 0 || jpix0 > BXM3) {return false;}
3181  if(ipix0 < 0 || ipix0 > BYM3) {return false;}
3182 
3183  xf = xhit + 0.5*zsize_*cota_current_ + lorxwidth_;
3184  yf = yhit + 0.5*zsize_*cotb_current_ + lorywidth_;
3185 
3186  jpixf = floor(xf/xsize_)+1;
3187  ipixf = floor(yf/ysize_)+1;
3188 
3189  if(jpixf < 0 || jpixf > BXM3) {return false;}
3190  if(ipixf < 0 || ipixf > BYM3) {return false;}
3191 
3192  // total charge length
3193 
3194  sf = std::sqrt((xf-x0)*(xf-x0) + (yf-y0)*(yf-y0));
3195  if((xf-x0) != 0.f) {slopey = (yf-y0)/(xf-x0);} else { slopey = 1.e10;}
3196  if((yf-y0) != 0.f) {slopex = (xf-x0)/(yf-y0);} else { slopex = 1.e10;}
3197 
3198  // use average charge in this direction
3199 
3200  qtotal = qavg_avg_;
3201 
3202  SimplePixel element;
3203  element.s = sf;
3204  element.x = xf;
3205  element.y = yf;
3206  element.i = ipixf;
3207  element.j = jpixf;
3208  element.btype = 0;
3209  list.push_back(element);
3210 
3211  // nx is the number of x interfaces crossed by the line charge
3212 
3213  nx = jpixf - jpix0;
3214  anx = abs(nx);
3215  if(anx > 0) {
3216  if(nx > 0) {
3217  for(j=jpix0; j<jpixf; ++j) {
3218  xi = xsize_*j;
3219  yi = slopey*(xi-x0) + y0;
3220  ipix = (int)(yi/ysize_)+1;
3221  si = std::sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
3222  element.s = si;
3223  element.x = xi;
3224  element.y = yi;
3225  element.i = ipix;
3226  element.j = j;
3227  element.btype = 1;
3228  list.push_back(element);
3229  }
3230  } else {
3231  for(j=jpix0; j>jpixf; --j) {
3232  xi = xsize_*(j-1);
3233  yi = slopey*(xi-x0) + y0;
3234  ipix = (int)(yi/ysize_)+1;
3235  si = std::sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
3236  element.s = si;
3237  element.x = xi;
3238  element.y = yi;
3239  element.i = ipix;
3240  element.j = j;
3241  element.btype = 1;
3242  list.push_back(element);
3243  }
3244  }
3245  }
3246 
3247  ny = ipixf - ipix0;
3248  any = abs(ny);
3249  if(any > 0) {
3250  if(ny > 0) {
3251  for(i=ipix0; i<ipixf; ++i) {
3252  yi = ysize_*i;
3253  xi = slopex*(yi-y0) + x0;
3254  jpix = (int)(xi/xsize_)+1;
3255  si = std::sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
3256  element.s = si;
3257  element.x = xi;
3258  element.y = yi;
3259  element.i = i;
3260  element.j = jpix;
3261  element.btype = 2;
3262  list.push_back(element);
3263  }
3264  } else {
3265  for(i=ipix0; i>ipixf; --i) {
3266  yi = ysize_*(i-1);
3267  xi = slopex*(yi-y0) + x0;
3268  jpix = (int)(xi/xsize_)+1;
3269  si = std::sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
3270  element.s = si;
3271  element.x = xi;
3272  element.y = yi;
3273  element.i = i;
3274  element.j = jpix;
3275  element.btype = 2;
3276  list.push_back(element);
3277  }
3278  }
3279  }
3280 
3281  imax = std::max(ipix0, ipixf);
3282  jmax = std::max(jpix0, jpixf);
3283 
3284  // Sort the list according to the distance from the initial point
3285 
3286  list.sort();
3287 
3288  // Look for double pixels and adjust the list appropriately
3289 
3290  for(i=1; i<imax; ++i) {
3291  if(ydouble[i-1]) {
3292  listIter = list.begin();
3293  if(ny > 0) {
3294  while(listIter != list.end()) {
3295  if(listIter->i == i && listIter->btype == 2) {
3296  listIter = list.erase(listIter);
3297  continue;
3298  }
3299  if(listIter->i > i) {
3300  --(listIter->i);
3301  }
3302  ++listIter;
3303  }
3304  } else {
3305  while(listIter != list.end()) {
3306  if(listIter->i == i+1 && listIter->btype == 2) {
3307  listIter = list.erase(listIter);
3308  continue;
3309  }
3310  if(listIter->i > i+1) {
3311  --(listIter->i);
3312  }
3313  ++listIter;
3314  }
3315  }
3316  }
3317  }
3318 
3319  for(j=1; j<jmax; ++j) {
3320  if(xdouble[j-1]) {
3321  listIter = list.begin();
3322  if(nx > 0) {
3323  while(listIter != list.end()) {
3324  if(listIter->j == j && listIter->btype == 1) {
3325  listIter = list.erase(listIter);
3326  continue;
3327  }
3328  if(listIter->j > j) {
3329  --(listIter->j);
3330  }
3331  ++listIter;
3332  }
3333  } else {
3334  while(listIter != list.end()) {
3335  if(listIter->j == j+1 && listIter->btype == 1) {
3336  listIter = list.erase(listIter);
3337  continue;
3338  }
3339  if(listIter->j > j+1) {
3340  --(listIter->j);
3341  }
3342  ++listIter;
3343  }
3344  }
3345  }
3346  }
3347 
3348  // The list now contains the path lengths of the line charge in each pixel from (x0,y0). Cacluate the lengths of the segments and the charge.
3349 
3350  s0 = 0.f;
3351  listIter = list.begin();
3352  listEnd = list.end();
3353  for( ;listIter != listEnd; ++listIter) {
3354  si = listIter->s;
3355  ds = si - s0;
3356  s0 = si;
3357  j = listIter->j;
3358  i = listIter->i;
3359  if(sf > 0.f) { qpix = qtotal*ds/sf;} else {qpix = qtotal;}
3360  template2d[j][i] += qpix;
3361  }
3362 
3363  return true;
3364 
3365 } // simpletemplate2D
3366 
3367 
3368 // ************************************************************************************************************
3373 // ************************************************************************************************************
3374 void SiPixelTemplate::vavilov_pars(double& mpv, double& sigma, double& kappa)
3375 
3376 {
3377  // Local variables
3378  int i;
3379  int ilow, ihigh, Ny;
3380  float yratio, cotb, cotalpha0, arg;
3381 
3382  // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta)
3383 
3384  cotalpha0 = thePixelTemp_[index_id_].enty[0].cotalpha;
3385  arg = cotb_current_*cotb_current_ + cota_current_*cota_current_ - cotalpha0*cotalpha0;
3386  if(arg < 0.f) arg = 0.f;
3387  cotb = std::sqrt(arg);
3388 
3389  // Copy the charge scaling factor to the private variable
3390 
3391  Ny = thePixelTemp_[index_id_].head.NTy;
3392 
3393 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3394  if(Ny < 2) {
3395  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny << std::endl;
3396  }
3397 #else
3398  assert(Ny > 1);
3399 #endif
3400 
3401  // next, loop over all y-angle entries
3402 
3403  ilow = 0;
3404  yratio = 0.f;
3405 
3406  if(cotb >= thePixelTemp_[index_id_].enty[Ny-1].cotbeta) {
3407 
3408  ilow = Ny-2;
3409  yratio = 1.f;
3410 
3411  } else {
3412 
3413  if(cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
3414 
3415  for (i=0; i<Ny-1; ++i) {
3416 
3417  if( thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i+1].cotbeta) {
3418 
3419  ilow = i;
3420  yratio = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta)/(thePixelTemp_[index_id_].enty[i+1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
3421  break;
3422  }
3423  }
3424  }
3425  }
3426 
3427  ihigh=ilow + 1;
3428 
3429  // Interpolate Vavilov parameters
3430 
3431  mpvvav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].mpvvav + yratio*thePixelTemp_[index_id_].enty[ihigh].mpvvav;
3432  sigmavav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sigmavav + yratio*thePixelTemp_[index_id_].enty[ihigh].sigmavav;
3433  kappavav_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].kappavav + yratio*thePixelTemp_[index_id_].enty[ihigh].kappavav;
3434 
3435  // Copy to parameter list
3436 
3437 
3438  mpv = (double)mpvvav_;
3439  sigma = (double)sigmavav_;
3440  kappa = (double)kappavav_;
3441 
3442  return;
3443 
3444 } // vavilov_pars
3445 
3446 // ************************************************************************************************************
3451 // ************************************************************************************************************
3452 void SiPixelTemplate::vavilov2_pars(double& mpv, double& sigma, double& kappa)
3453 
3454 {
3455  // Local variables
3456  int i;
3457  int ilow, ihigh, Ny;
3458  float yratio, cotb, cotalpha0, arg;
3459 
3460  // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta)
3461 
3462  cotalpha0 = thePixelTemp_[index_id_].enty[0].cotalpha;
3463  arg = cotb_current_*cotb_current_ + cota_current_*cota_current_ - cotalpha0*cotalpha0;
3464  if(arg < 0.f) arg = 0.f;
3465  cotb = std::sqrt(arg);
3466 
3467  // Copy the charge scaling factor to the private variable
3468 
3469  Ny = thePixelTemp_[index_id_].head.NTy;
3470 
3471 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
3472  if(Ny < 2) {
3473  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny << std::endl;
3474  }
3475 #else
3476  assert(Ny > 1);
3477 #endif
3478 
3479  // next, loop over all y-angle entries
3480 
3481  ilow = 0;
3482  yratio = 0.f;
3483 
3484  if(cotb >= thePixelTemp_[index_id_].enty[Ny-1].cotbeta) {
3485 
3486  ilow = Ny-2;
3487  yratio = 1.f;
3488 
3489  } else {
3490 
3491  if(cotb >= thePixelTemp_[index_id_].enty[0].cotbeta) {
3492 
3493  for (i=0; i<Ny-1; ++i) {
3494 
3495  if( thePixelTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < thePixelTemp_[index_id_].enty[i+1].cotbeta) {
3496 
3497  ilow = i;
3498  yratio = (cotb - thePixelTemp_[index_id_].enty[i].cotbeta)/(thePixelTemp_[index_id_].enty[i+1].cotbeta - thePixelTemp_[index_id_].enty[i].cotbeta);
3499  break;
3500  }
3501  }
3502  }
3503  }
3504 
3505  ihigh=ilow + 1;
3506 
3507  // Interpolate Vavilov parameters
3508 
3509  mpvvav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].mpvvav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].mpvvav2;
3510  sigmavav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].sigmavav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].sigmavav2;
3511  kappavav2_ = (1.f - yratio)*thePixelTemp_[index_id_].enty[ilow].kappavav2 + yratio*thePixelTemp_[index_id_].enty[ihigh].kappavav2;
3512 
3513  // Copy to parameter list
3514 
3515  mpv = (double)mpvvav2_;
3516  sigma = (double)sigmavav2_;
3517  kappa = (double)kappavav2_;
3518 
3519  return;
3520 
3521 } // vavilov2_pars
3522 
3523 
SiPixelTemplateEntry::qavg_avg
float qavg_avg
average cluster charge of clusters that are less than qavg (normalize 2-D simple templates) ...
Definition: SiPixelTemplate.h:179
SiPixelTemplateEntry::runnum
int runnum
< Basic template entry corresponding to a single set of track angles
Definition: SiPixelTemplate.h:107
SiPixelTemplateEntry::xflpar
float xflpar[4][6]
Aqfl-parameterized x-correction in 4 charge bins.
Definition: SiPixelTemplate.h:153
SiPixelTemplateEntry::xrms
float xrms[4]
average x-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:150
SiPixelTemplate::vavilov2_pars
void vavilov2_pars(double &mpv, double &sigma, double &kappa)
Definition: SiPixelTemplate.cc:3452
BXSIZE
#define BXSIZE
Definition: SiPixelTemplateDefs.h:33
LOGERROR
#define LOGERROR(x)
Definition: SiPixelTemplate.cc:112
SiPixelTemplateEntry::clslenx
float clslenx
cluster x-length in pixels at signal height sxmax/2
Definition: SiPixelTemplate.h:133
SiPixelTemplateEntry::qavg
float qavg
average cluster charge for this set of track angles (now includes threshold effects) ...
Definition: SiPixelTemplate.h:113
mps_fire.i
i
Definition: mps_fire.py:338
SiPixelTemplateEntry::syone
float syone
rms for one pixel y-clusters
Definition: SiPixelTemplate.h:117
SiPixelTemplateEntry::chi2yavgone
float chi2yavgone
average y chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:158
SiPixelTemplate::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
Definition: SiPixelTemplate.cc:1001
SiPixelTemplate::cytemp
int cytemp()
Return central pixel of y template pixels above readout threshold.
Definition: SiPixelTemplate.cc:2058
SiPixelTemplateEntry::dyone
float dyone
mean offset/correction for one pixel y-clusters
Definition: SiPixelTemplate.h:116
BYSIZE
#define BYSIZE
Definition: SiPixelTemplateDefs.h:25
SiPixelTemplateEntry::sigmavav
float sigmavav
"sigma" scale fctor for Vavilov distribution
Definition: SiPixelTemplate.h:135
SiPixelTemplateEntry::fracxtwo
float fracxtwo
fraction of double pixel sample with xsize = 1
Definition: SiPixelTemplate.h:178
SiPixelTemplateEntry::yavggen
float yavggen[4]
generic algorithm: average y-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:127
SiPixelTemplateEntry::sxmax
float sxmax
average pixel signal for x-projection of cluster
Definition: SiPixelTemplate.h:118
SimplePixel::btype
int btype
type of boundary (0=end, 1 = x-boundary, 2 = y-boundary)
Definition: SimplePixel.h:22
SiPixelTemplate::xsigma2
void xsigma2(int fxpix, int lxpix, float sxthr, float xsum[13+4], float xsig2[13+4])
Definition: SiPixelTemplate.cc:1691
SiPixelTemplateEntry::chi2yminc2m
float chi2yminc2m[4]
1st pass chi2 min search: minimum of y chi^2 in 4 charge bins (merged clusters)
Definition: SiPixelTemplate.h:165
SiPixelTemplateEntry::xavgc2m
float xavgc2m[4]
1st pass chi2 min search: average x-bias of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:166
TXSIZE
#define TXSIZE
Definition: SiPixelTemplateDefs.h:30
any
bool any(const std::vector< T > &v, const T &what)
Definition: ECalSD.cc:37
SiPixelTemplateEntry::chi2xavgone
float chi2xavgone
average x chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:160
SiPixelTemplateEntry::fracytwo
float fracytwo
fraction of double pixel sample with ysize = 1
Definition: SiPixelTemplate.h:177
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
SiPixelTemplateEntry::yrms
float yrms[4]
average y-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:145
SiPixelTemplateHeader::NTxx
int NTxx
number of Template x-entries in each slice
Definition: SiPixelTemplate.h:191
BXM3
#define BXM3
Definition: SiPixelTemplateDefs.h:37
BXM1
#define BXM1
Definition: SiPixelTemplateDefs.h:35
SiPixelTemplateEntry::mpvvav2
float mpvvav2
most probable charge in Vavilov distribution for 2 merged clusters (not actually for larger kappa) ...
Definition: SiPixelTemplate.h:137
SiPixelTemplateEntry::ygx0gen
float ygx0gen[4]
generic algorithm: average y0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:170
SiPixelTemplateEntry::xavg
float xavg[4]
average x-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:149
SiPixelTemplate::simpletemplate2D
bool simpletemplate2D(float xhitp, float yhitp, std::vector< bool > &ydouble, std::vector< bool > &xdouble, float template2d[13+2][21+2])
Make simple 2-D templates from track angles set in interpolate and hit position.
Definition: SiPixelTemplate.cc:3160
SiPixelTemplateHeader::lorybias
float lorybias
estimate of y-lorentz bias
Definition: SiPixelTemplate.h:196
SimplePixel::y
float y
y coordinate of boundary intersection
Definition: SimplePixel.h:19
SiPixelTemplateHeader::fluence
float fluence
radiation fluence in n_eq/cm^2
Definition: SiPixelTemplate.h:200
SiPixelTemplateHeader::ss50
float ss50
1/2 of the single hit dcol threshold in electrons
Definition: SiPixelTemplate.h:202
SiPixelTemplateEntry::dytwo
float dytwo
mean offset/correction for one double-pixel y-clusters
Definition: SiPixelTemplate.h:121
SiPixelTemplateHeader::ysize
float ysize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate.h:208
SiPixelTemplateEntry::dxone
float dxone
mean offset/correction for one pixel x-clusters
Definition: SiPixelTemplate.h:119
SiPixelTemplateEntry::pixmax
float pixmax
maximum charge for individual pixels in cluster
Definition: SiPixelTemplate.h:114
SiPixelTemplateHeader::qscale
float qscale
Charge scaling to match cmssw and pixelav.
Definition: SiPixelTemplate.h:193
SiPixelTemplateEntry::chi2yavgc2m
float chi2yavgc2m[4]
1st pass chi2 min search: average y chi^2 in 4 charge bins (merged clusters)
Definition: SiPixelTemplate.h:164
SiPixelTemplateEntry::qmin
float qmin
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
Definition: SiPixelTemplate.h:125
SiPixelTemplateEntry::chi2xavgc2m
float chi2xavgc2m[4]
1st pass chi2 min search: average x chi^2 in 4 charge bins (merged clusters)
Definition: SiPixelTemplate.h:168
SiPixelTemplate::xflcorr
float xflcorr(int binq, float qflx)
Definition: SiPixelTemplate.cc:1854
SiPixelTemplateEntry::xgsiggen
float xgsiggen[4]
generic algorithm: average sigma_x from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:173
protons_cff.xi
xi
Definition: protons_cff.py:31
SiPixelTemplateEntry::chi2xminc2m
float chi2xminc2m[4]
1st pass chi2 min search: minimum of x chi^2 in 4 charge bins (merged clusters)
Definition: SiPixelTemplate.h:169
SiPixelTemplateHeader::ID
int ID
< template header structure
Definition: SiPixelTemplate.h:188
edm::isFinite
constexpr bool isFinite(T x)
SiPixelTemplateHeader::Bfield
float Bfield
Bfield in Tesla.
Definition: SiPixelTemplate.h:205
SiPixelTemplateEntry::cotalpha
float cotalpha
cot(alpha) is proportional to cluster length in x and is basis of interpolation
Definition: SiPixelTemplate.h:109
SiPixelTemplateHeader::temperature
float temperature
detector temperature in deg K
Definition: SiPixelTemplate.h:199
funct::arg
A arg
Definition: Factorize.h:38
SiPixelTemplateEntry::ytemp
float ytemp[9][21]
templates for y-reconstruction (binned over 1 central pixel)
Definition: SiPixelTemplate.h:141
SiPixelTemplateEntry::yrmsc2m
float yrmsc2m[4]
1st pass chi2 min search: average y-rms of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:163
SiPixelTemplateHeader::lorywidth
float lorywidth
estimate of y-lorentz width for optimal resolution
Definition: SiPixelTemplate.h:194
BYM1
#define BYM1
Definition: SiPixelTemplateDefs.h:27
BXM2
#define BXM2
Definition: SiPixelTemplateDefs.h:36
SiPixelTemplateHeader::fbin
float fbin[3]
The QBin definitions in Q_clus/Q_avg.
Definition: SiPixelTemplate.h:206
SiPixelTemplate::ytemp
void ytemp(int fybin, int lybin, float ytemplate[41][21+4])
Definition: SiPixelTemplate.cc:1906
SiPixelTemplate::cxtemp
int cxtemp()
Return central pixel of x-template pixels above readout threshold.
Definition: SiPixelTemplate.cc:2107
pfDeepBoostedJetPreprocessParams_cfi.lower_bound
lower_bound
Definition: pfDeepBoostedJetPreprocessParams_cfi.py:15
SiPixelTemplateHeader::templ_version
int templ_version
Version number of the template to ensure code compatibility.
Definition: SiPixelTemplate.h:204
SiPixelTemplateEntry::xpar
float xpar[2][5]
projected x-pixel uncertainty parameterization
Definition: SiPixelTemplate.h:142
SiPixelTemplateEntry::ygx0
float ygx0[4]
average y0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:146
SiPixelTemplateEntry::xgx0gen
float xgx0gen[4]
generic algorithm: average x0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:172
SiPixelTemplateEntry::symax
float symax
average pixel signal for y-projection of cluster
Definition: SiPixelTemplate.h:115
SiPixelTemplateEntry::kappavav2
float kappavav2
kappa parameter for Vavilov distribution for 2 merged clusters
Definition: SiPixelTemplate.h:139
SiPixelTemplateEntry::xgsig
float xgsig[4]
average sigma_x from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:152
SiPixelTemplateEntry::cotbeta
float cotbeta
cot(beta) is proportional to cluster length in y and is basis of interpolation
Definition: SiPixelTemplate.h:111
SiPixelTemplateEntry::beta
float beta
beta track angle (defined in CMS CMS IN 2004/014)
Definition: SiPixelTemplate.h:110
SiPixelTemplate::vavilov_pars
void vavilov_pars(double &mpv, double &sigma, double &kappa)
Definition: SiPixelTemplate.cc:3374
triggerObjects_cff.id
id
Definition: triggerObjects_cff.py:29
ENDL
#define ENDL
Definition: SiPixelTemplate.cc:115
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
BHX
#define BHX
Definition: SiPixelTemplateDefs.h:34
SiPixelTemplate::yflcorr
float yflcorr(int binq, float qfly)
Definition: SiPixelTemplate.cc:1802
SimplePixel::s
float s
distance from track entry
Definition: SimplePixel.h:17
SiPixelTemplateEntry::xrmsc2m
float xrmsc2m[4]
1st pass chi2 min search: average x-rms of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:167
FrontierConditions_GlobalTag_cff.file
file
Definition: FrontierConditions_GlobalTag_cff.py:13
SiPixelTemplate::ytemp3d
void ytemp3d(int j, int k, std::vector< float > &ytemplate)
Definition: SiPixelTemplate.cc:2235
SiPixelTemplate::temperrors
void temperrors(int id, float cotalpha, float cotbeta, int qBin, float &sigmay, float &sigmax, float &sy1, float &sy2, float &sx1, float &sx2)
Definition: SiPixelTemplate.cc:2777
SiPixelTemplateHeader::lorxwidth
float lorxwidth
estimate of x-lorentz width for optimal resolution
Definition: SiPixelTemplate.h:195
SiPixelTemplateEntry::xgx0
float xgx0[4]
average x0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:151
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
SiPixelTemplate::ysigma2
void ysigma2(int fypix, int lypix, float sythr, float ysum[21+4], float ysig2[21+4])
SiPixelTemplateHeader::title
char title[80]
template title
Definition: SiPixelTemplate.h:203
BYM2
#define BYM2
Definition: SiPixelTemplateDefs.h:28
SimplePixel::j
int j
y index of traversed pixel
Definition: SimplePixel.h:21
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
SiPixelTemplateEntry::ygsiggen
float ygsiggen[4]
generic algorithm: average sigma_y from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:171
SiPixelTemplateEntry::fracxone
float fracxone
fraction of sample with xsize = 1
Definition: SiPixelTemplate.h:176
SiPixelTemplateEntry::mpvvav
float mpvvav
most probable charge in Vavilov distribution (not actually for larger kappa)
Definition: SiPixelTemplate.h:134
BYM3
#define BYM3
Definition: SiPixelTemplateDefs.h:29
SiPixelTemplateEntry::sxtwo
float sxtwo
rms for one double-pixel x-clusters
Definition: SiPixelTemplate.h:124
SiPixelTemplateEntry::r_qMeas_qTrue
float r_qMeas_qTrue
ratio of measured to true cluster charge
Definition: SiPixelTemplate.h:180
SiPixelTemplateHeader::Dtype
int Dtype
detector type (0=BPix, 1=FPix)
Definition: SiPixelTemplate.h:192
SiPixelTemplateEntry::chi2yavg
float chi2yavg[4]
average y chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:154
SiPixelTemplateEntry::alpha
float alpha
alpha track angle (defined in CMS CMS IN 2004/014)
Definition: SiPixelTemplate.h:108
SiPixelTemplateEntry::xtemp
float xtemp[9][13]
templates for x-reconstruction (binned over 1 central pixel)
Definition: SiPixelTemplate.h:143
SiPixelTemplateStore::head
SiPixelTemplateHeader head
< template storage structure
Definition: SiPixelTemplate.h:215
SiPixelTemplateEntry::chi2xminone
float chi2xminone
minimum of x chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:161
SiPixelTemplateEntry::ypar
float ypar[2][5]
projected y-pixel uncertainty parameterization
Definition: SiPixelTemplate.h:140
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
SiPixelTemplate::xtemp3d
void xtemp3d(int j, int k, std::vector< float > &xtemplate)
Definition: SiPixelTemplate.cc:2338
SiPixelTemplateHeader::s50
float s50
1/2 of the multihit dcol threshold in electrons
Definition: SiPixelTemplate.h:201
TYSIZE
#define TYSIZE
Definition: SiPixelTemplateDefs.h:21
SiPixelTemplate::xtemp
void xtemp(int fxbin, int lxbin, float xtemplate[41][13+4])
Definition: SiPixelTemplate.cc:1984
SiPixelTemplateHeader::zsize
float zsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate.h:209
SiPixelTemplate::qbin_dist
void qbin_dist(int id, float cotalpha, float cotbeta, float qbin_frac[4], float &ny1_frac, float &ny2_frac, float &nx1_frac, float &nx2_frac)
Definition: SiPixelTemplate.cc:2976
SiPixelTemplateHeader::Vbias
float Vbias
detector bias potential in Volts
Definition: SiPixelTemplate.h:198
SiPixelTemplateEntry::qbfrac
float qbfrac[3]
fraction of sample in qbin = 0-2 (>=3 is the complement)
Definition: SiPixelTemplate.h:174
SiPixelTemplateHeader::xsize
float xsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate.h:207
SiPixelTemplate::ytemp3d_int
void ytemp3d_int(int nypix, int &nybins)
Definition: SiPixelTemplate.cc:2157
SiPixelTemplateEntry::fracyone
float fracyone
fraction of sample with ysize = 1
Definition: SiPixelTemplate.h:175
SiPixelTemplateEntry::qmin2
float qmin2
tighter minimum cluster charge for valid hit (keeps 99.8% of simulated hits)
Definition: SiPixelTemplate.h:126
SiPixelTemplateEntry::yavg
float yavg[4]
average y-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:144
SiPixelTemplate::xtemp3d_int
void xtemp3d_int(int nxpix, int &nxbins)
Definition: SiPixelTemplate.cc:2259
BHY
#define BHY
Definition: SiPixelTemplateDefs.h:26
SiPixelTemplate::pushfile
static bool pushfile(int filenum, std::vector< SiPixelTemplateStore > &pixelTemp, std::string dir="CalibTracker/SiPixelESProducers/data/")
Definition: SiPixelTemplate.cc:132
SiPixelTemplateEntry::chi2ymin
float chi2ymin[4]
minimum of y chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:155
SiPixelTemplateEntry::yavgc2m
float yavgc2m[4]
1st pass chi2 min search: average y-bias of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:162
SimplePixel::x
float x
x coordinate of boundary intersection
Definition: SimplePixel.h:18
SiPixelTemplateHeader::NTy
int NTy
number of Template y entries
Definition: SiPixelTemplate.h:189
SiPixelTemplateEntry::dxtwo
float dxtwo
mean offset/correction for one double-pixel x-clusters
Definition: SiPixelTemplate.h:123
SiPixelTemplateHeader::lorxbias
float lorxbias
estimate of x-lorentz bias
Definition: SiPixelTemplate.h:197
edm
HLT enums.
Definition: AlignableModifier.h:17
LOGINFO
#define LOGINFO(x)
Definition: SiPixelTemplate.cc:113
SimplePixel::i
int i
x index of traversed pixel
Definition: SimplePixel.h:20
SiPixelTemplateHeader::NTyx
int NTyx
number of Template y-slices of x entries
Definition: SiPixelTemplate.h:190
SiPixelTemplateEntry::yflpar
float yflpar[4][6]
Aqfl-parameterized y-correction in 4 charge bins.
Definition: SiPixelTemplate.h:148
edm::FileInPath::fullPath
std::string fullPath() const
Definition: FileInPath.cc:163
SiPixelTemplateEntry::sxone
float sxone
rms for one pixel x-clusters
Definition: SiPixelTemplate.h:120
ApeEstimator_cff.qBin
qBin
Definition: ApeEstimator_cff.py:34
SiPixelTemplateEntry::chi2xmin
float chi2xmin[4]
minimum of x chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:157
kappa
static const G4double kappa
Definition: UrbanMscModel93.cc:35
gather_cfg.cout
cout
Definition: gather_cfg.py:144
SiPixelTemplateEntry::chi2yminone
float chi2yminone
minimum of y chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:159
SiPixelTemplate::qbin
int qbin(int id, float cotalpha, float cotbeta, float locBz, float locBx, float qclus, float &pixmx, float &sigmay, float &deltay, float &sigmax, float &deltax, float &sy1, float &dy1, float &sy2, float &dy2, float &sx1, float &dx1, float &sx2, float &dx2)
Definition: SiPixelTemplate.cc:2386
dir
dbl *** dir
Definition: mlp_gen.cc:35
SiPixelTemplateEntry::costrk
float costrk[3]
direction cosines of tracks used to generate this entry
Definition: SiPixelTemplate.h:112
SiPixelTemplate::postInit
static void postInit(std::vector< SiPixelTemplateStore > &thePixelTemp_)
Definition: SiPixelTemplate.cc:960
SiPixelTemplateStore::entx
SiPixelTemplateEntry entx[5][29]
29 Barrel x templates spanning cluster lengths from -6px (-1.125Rad) to +6px (+1.125Rad) in each of 5...
Definition: SiPixelTemplate.h:221
SiPixelTemplateEntry::ygsig
float ygsig[4]
average sigma_y from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:147
SiPixelTemplateStore::enty
SiPixelTemplateEntry enty[60]
60 Barrel y templates spanning cluster lengths from 0px to +18px [28 entries for fpix] ...
Definition: SiPixelTemplate.h:220
SiPixelTemplateEntry::sigmavav2
float sigmavav2
"sigma" scale fctor for Vavilov distribution for 2 merged clusters
Definition: SiPixelTemplate.h:138
SiPixelTemplateEntry::sytwo
float sytwo
rms for one double-pixel y-clusters
Definition: SiPixelTemplate.h:122
SiPixelTemplateEntry::xrmsgen
float xrmsgen[4]
generic algorithm: average x-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:130
SiPixelTemplateEntry::chi2xavg
float chi2xavg[4]
average x chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:156
list
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger list("!*","!HLTx*"if it matches 2 triggers or more) will accept the event if all the matching triggers are FAIL.It will reject the event if any of the triggers are PASS or EXCEPTION(this matches the behavior of"!*"before the partial wildcard feature was incorporated).Triggers which are in the READY state are completely ignored.(READY should never be returned since the trigger paths have been run
SiPixelTemplateEntry::kappavav
float kappavav
kappa parameter for Vavilov distribution
Definition: SiPixelTemplate.h:136
SiPixelTemplateEntry::xavggen
float xavggen[4]
generic algorithm: average x-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:129
SiPixelTemplateEntry::yrmsgen
float yrmsgen[4]
generic algorithm: average y-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:128
SiPixelTemplateEntry::clsleny
float clsleny
cluster y-length in pixels at signal height symax/2
Definition: SiPixelTemplate.h:132
Generated for CMSSW Reference Manual by doxygen 1.8.11