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