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