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