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