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