CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiPixelTemplate.cc
Go to the documentation of this file.
1 //
2 // SiPixelTemplate.cc Version 8.13
3 //
4 // Add goodness-of-fit info and spare entries to templates, version number in template header, more error checking
5 // Add correction for (Q_F-Q_L)/(Q_F+Q_L) bias
6 // Add cot(beta) reflection to reduce y-entries and more sophisticated x-interpolation
7 // Fix small index searching bug in interpolate method
8 // Change interpolation indexing to avoid complier complaining about possible un-initialized variables
9 // Replace containers with static arrays in calls to ysigma2 and xsigma2
10 // Add external threshold to calls to ysigma2 and xsigma2, fix parameter signal max for xsigma2
11 // Return to 5 pixel spanning but adjust boundaries to use only when needed
12 // Implement improved (faster) chi2min search that depends on pixel types
13 // Fill template arrays in single calls to this object
14 // Add qmin to the template
15 // Add qscale to match charge scales
16 // Small improvement to x-chisquare interpolation
17 // Enlarge SiPixelTemplateStore to accommodate larger templates and increased alpha acceptance (reduce PT threshold to ~200 MeV)
18 // Change output streams to conform to CMSSW info and error logging.
19 // Store x and y cluster sizes in fractional pixels to facilitate cluster splitting
20 // Add methods to return 3-d templates needed for cluster splitting
21 // Keep interpolated central 9 template bins in private storage and expand/shift in the getter functions (faster for speed=2/3) and easier to build 3d templates
22 // Store error and bias information for the simple chi^2 min position analysis (no interpolation or Q_{FB} corrections) to use in cluster splitting
23 // To save time, the gaussian centers and sigma are not interpolated right now (they aren't currently used). They can be restored by un-commenting lines in the interpolate method.
24 // Add a new method to calculate qbin for input cotbeta and cluster charge. To be used for error estimation of merged clusters in PixelCPEGeneric.
25 // Improve the charge estimation for larger cot(alpha) tracks
26 // Change interpolate method to return false boolean if track angles are outside of range
27 // Add template info and method for truncation information
28 // Change to allow template sizes to be changed at compile time
29 // Fix bug in track angle checking
30 // Accommodate Dave's new DB pushfile which overloads the old method (file input)
31 // Add CPEGeneric error information and expand qbin method to access useful info for PixelCPEGeneric
32 // Fix large cot(alpha) bug in qmin interpolation
33 // Add second qmin to allow a qbin=5 state
34 // Use interpolated chi^2 info for one-pixel clusters
35 // Fix DB pushfile version number checking bug.
36 // Remove assert from qbin method
37 // Replace asserts with exceptions in CMSSW
38 // Change calling sequence to interpolate method to handle cot(beta)<0 for FPix cosmics
39 // Add getter for pixelav Lorentz width estimates to qbin method
40 // Add check on template size to interpolate and qbin methods
41 // Add qbin population information, charge distribution information
42 //
43 //
44 // V7.00 - Decouple BPix and FPix information into separate templates
45 // Add methods to facilitate improved cluster splitting
46 // Fix small charge scaling bug (affects FPix only)
47 // Change y-slice used for the x-template to be closer to the actual cotalpha-cotbeta point
48 // (there is some weak breakdown of x-y factorization in the FPix after irradiation)
49 //
50 //
51 // V8.00 - Add method to calculate a simple 2D template
52 // Reorganize the interpolate method to extract header info only once per ID
53 // V8.01 - Improve simple template normalization
54 // V8.05 - Change qbin normalization to work better after irradiation
55 // V8.10 - Add Vavilov distribution interpolation
56 // V8.11 - Renormalize the x-templates for Guofan's cluster size calculation
57 // V8.12 - Technical fix to qavg issue.
58 // V8.13 - Fix qbin and fastsim interpolaters to avoid changing class variables
59 
60 // Created by Morris Swartz on 10/27/06.
61 // Copyright 2006 __TheJohnsHopkinsUniversity__. All rights reserved.
62 //
63 //
64 
65 //#include <stdlib.h>
66 //#include <stdio.h>
67 #include <cmath>
68 #include <algorithm>
69 #include <vector>
70 //#include "boost/multi_array.hpp"
71 #include <iostream>
72 #include <iomanip>
73 #include <sstream>
74 #include <fstream>
75 #include <list>
76 
77 
78 
79 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
80 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplate.h"
81 #include "RecoLocalTracker/SiPixelRecHits/interface/SimplePixel.h"
82 #include "FWCore/ParameterSet/interface/FileInPath.h"
83 #include "FWCore/MessageLogger/interface/MessageLogger.h"
84 #define LOGERROR(x) LogError(x)
85 #define LOGINFO(x) LogInfo(x)
86 #define ENDL " "
87 #include "FWCore/Utilities/interface/Exception.h"
88 using namespace edm;
89 #else
90 #include "SiPixelTemplate.h"
91 #include "SimplePixel.h"
92 #define LOGERROR(x) std::cout << x << ": "
93 #define LOGINFO(x) std::cout << x << ": "
94 #define ENDL std::endl
95 #endif
96 
97 //****************************************************************
102 //****************************************************************
103 bool SiPixelTemplate::pushfile(int filenum)
104 {
105  // Add template stored in external file numbered filenum to theTemplateStore
106 
107  // Local variables
108  int i, j, k, l;
109  float qavg_avg;
110  const char *tempfile;
111  // char title[80]; remove this
112  char c;
113  const int code_version={16};
114 
115 
116 
117  // Create a filename for this run
118 
119  std::ostringstream tout;
120 
121  // Create different path in CMSSW than standalone
122 
123 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
124  tout << "CalibTracker/SiPixelESProducers/data/template_summary_zp"
125  << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
126  std::string tempf = tout.str();
127  edm::FileInPath file( tempf.c_str() );
128  tempfile = (file.fullPath()).c_str();
129 #else
130  tout << "template_summary_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
131  std::string tempf = tout.str();
132  tempfile = tempf.c_str();
133 #endif
134 
135  // open the template file
136 
137  std::ifstream in_file(tempfile, std::ios::in);
138 
139  if(in_file.is_open()) {
140 
141  // Create a local template storage entry
142 
143  SiPixelTemplateStore theCurrentTemp;
144 
145  // Read-in a header string first and print it
146 
147  for (i=0; (c=in_file.get()) != '\n'; ++i) {
148  if(i < 79) {theCurrentTemp.head.title[i] = c;}
149  }
150  if(i > 78) {i=78;}
151  theCurrentTemp.head.title[i+1] ='\0';
152  LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
153 
154  // next, the header information
155 
156  in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
157  >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
158  >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
159 
160  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file, no template load" << ENDL; return false;}
161 
162  LOGINFO("SiPixelTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
163  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
164  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
165  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
166  << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
167  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
168 
169  if(theCurrentTemp.head.templ_version < code_version) {LOGERROR("SiPixelTemplate") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
170 
171  // next, loop over all y-angle entries
172 
173  for (i=0; i < theCurrentTemp.head.NTy; ++i) {
174 
175  in_file >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0]
176  >> theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2];
177 
178  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
179 
180  // Calculate the alpha, beta, and cot(beta) for this entry
181 
182  theCurrentTemp.enty[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
183 
184  theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0]/theCurrentTemp.enty[i].costrk[2];
185 
186  theCurrentTemp.enty[i].beta = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
187 
188  theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1]/theCurrentTemp.enty[i].costrk[2];
189 
190  in_file >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].symax >> theCurrentTemp.enty[i].dyone
191  >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].sxmax >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
192 
193  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
194 
195  in_file >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo
196  >> theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clsleny >> theCurrentTemp.enty[i].clslenx;
197 
198  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 3, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
199 
200  for (j=0; j<2; ++j) {
201 
202  in_file >> theCurrentTemp.enty[i].ypar[j][0] >> theCurrentTemp.enty[i].ypar[j][1]
203  >> theCurrentTemp.enty[i].ypar[j][2] >> theCurrentTemp.enty[i].ypar[j][3] >> theCurrentTemp.enty[i].ypar[j][4];
204 
205  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 4, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
206 
207  }
208 
209  for (j=0; j<9; ++j) {
210 
211  for (k=0; k<TYSIZE; ++k) {in_file >> theCurrentTemp.enty[i].ytemp[j][k];}
212 
213  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 5, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
214  }
215 
216  for (j=0; j<2; ++j) {
217 
218  in_file >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1]
219  >> theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
220 
221  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 6, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
222 
223  }
224 
225  qavg_avg = 0.;
226  for (j=0; j<9; ++j) {
227 
228  for (k=0; k<TXSIZE; ++k) {in_file >> theCurrentTemp.enty[i].xtemp[j][k]; qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];}
229 
230  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 7, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
231  }
232  theCurrentTemp.enty[i].qavg_avg = qavg_avg/9.;
233 
234  for (j=0; j<4; ++j) {
235 
236  in_file >> theCurrentTemp.enty[i].yavg[j] >> theCurrentTemp.enty[i].yrms[j] >> theCurrentTemp.enty[i].ygx0[j] >> theCurrentTemp.enty[i].ygsig[j];
237 
238  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 8, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
239  }
240 
241  for (j=0; j<4; ++j) {
242 
243  in_file >> theCurrentTemp.enty[i].yflpar[j][0] >> theCurrentTemp.enty[i].yflpar[j][1] >> theCurrentTemp.enty[i].yflpar[j][2]
244  >> theCurrentTemp.enty[i].yflpar[j][3] >> theCurrentTemp.enty[i].yflpar[j][4] >> theCurrentTemp.enty[i].yflpar[j][5];
245 
246  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 9, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
247  }
248 
249  for (j=0; j<4; ++j) {
250 
251  in_file >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >> theCurrentTemp.enty[i].xgsig[j];
252 
253  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 10, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
254  }
255 
256  for (j=0; j<4; ++j) {
257 
258  in_file >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >> theCurrentTemp.enty[i].xflpar[j][2]
259  >> theCurrentTemp.enty[i].xflpar[j][3] >> theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
260 
261  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 11, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
262  }
263 
264  for (j=0; j<4; ++j) {
265 
266  in_file >> theCurrentTemp.enty[i].chi2yavg[j] >> theCurrentTemp.enty[i].chi2ymin[j] >> theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j];
267 
268  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 12, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
269  }
270 
271  for (j=0; j<4; ++j) {
272 
273  in_file >> theCurrentTemp.enty[i].yavgc2m[j] >> theCurrentTemp.enty[i].yrmsc2m[j] >> theCurrentTemp.enty[i].ygx0c2m[j] >> theCurrentTemp.enty[i].ygsigc2m[j];
274 
275  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 13, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
276  }
277 
278  for (j=0; j<4; ++j) {
279 
280  in_file >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >> theCurrentTemp.enty[i].xgx0c2m[j] >> theCurrentTemp.enty[i].xgsigc2m[j];
281 
282  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
283  }
284 
285  for (j=0; j<4; ++j) {
286 
287  in_file >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >> theCurrentTemp.enty[i].ygx0gen[j] >> theCurrentTemp.enty[i].ygsiggen[j];
288 
289  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14a, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
290  }
291 
292  for (j=0; j<4; ++j) {
293 
294  in_file >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >> theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
295 
296  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14b, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
297  }
298 
299  in_file >> theCurrentTemp.enty[i].chi2yavgone >> theCurrentTemp.enty[i].chi2yminone >> theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2
300  >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav >> theCurrentTemp.enty[i].qavg_spare >> theCurrentTemp.enty[i].spare[0];
301 
302  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 15, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
303 
304  in_file >> theCurrentTemp.enty[i].spare[1] >> theCurrentTemp.enty[i].spare[2] >> theCurrentTemp.enty[i].spare[3] >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1]
305  >> theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracyone >> theCurrentTemp.enty[i].fracxone >> theCurrentTemp.enty[i].fracytwo >> theCurrentTemp.enty[i].fracxtwo;
306  // theCurrentTemp.enty[i].qbfrac[3] = 1. - theCurrentTemp.enty[i].qbfrac[0] - theCurrentTemp.enty[i].qbfrac[1] - theCurrentTemp.enty[i].qbfrac[2];
307 
308  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 16, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
309 
310  }
311 
312  // next, loop over all barrel x-angle entries
313 
314  for (k=0; k < theCurrentTemp.head.NTyx; ++k) {
315 
316  for (i=0; i < theCurrentTemp.head.NTxx; ++i) {
317 
318  in_file >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0]
319  >> theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2];
320 
321  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 17, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
322 
323  // Calculate the alpha, beta, and cot(beta) for this entry
324 
325  theCurrentTemp.entx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
326 
327  theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0]/theCurrentTemp.entx[k][i].costrk[2];
328 
329  theCurrentTemp.entx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
330 
331  theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1]/theCurrentTemp.entx[k][i].costrk[2];
332 
333  in_file >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >> theCurrentTemp.entx[k][i].symax >> theCurrentTemp.entx[k][i].dyone
334  >> theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
335 
336  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 18, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
337 
338  in_file >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >> theCurrentTemp.entx[k][i].dxtwo
339  >> theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clsleny >> theCurrentTemp.entx[k][i].clslenx;
340  // >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
341 
342  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 19, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
343 
344  for (j=0; j<2; ++j) {
345 
346  in_file >> theCurrentTemp.entx[k][i].ypar[j][0] >> theCurrentTemp.entx[k][i].ypar[j][1]
347  >> theCurrentTemp.entx[k][i].ypar[j][2] >> theCurrentTemp.entx[k][i].ypar[j][3] >> theCurrentTemp.entx[k][i].ypar[j][4];
348 
349  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 20, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
350  }
351 
352  for (j=0; j<9; ++j) {
353 
354  for (l=0; l<TYSIZE; ++l) {in_file >> theCurrentTemp.entx[k][i].ytemp[j][l];}
355 
356  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 21, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
357  }
358 
359  for (j=0; j<2; ++j) {
360 
361  in_file >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1]
362  >> theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >> theCurrentTemp.entx[k][i].xpar[j][4];
363 
364 
365  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 22, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
366  }
367 
368  qavg_avg = 0.;
369  for (j=0; j<9; ++j) {
370 
371  for (l=0; l<TXSIZE; ++l) {in_file >> theCurrentTemp.entx[k][i].xtemp[j][l]; qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];}
372 
373  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 23, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
374  }
375  theCurrentTemp.entx[k][i].qavg_avg = qavg_avg/9.;
376 
377  for (j=0; j<4; ++j) {
378 
379  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];
380 
381  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 24, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
382  }
383 
384  for (j=0; j<4; ++j) {
385 
386  in_file >> theCurrentTemp.entx[k][i].yflpar[j][0] >> theCurrentTemp.entx[k][i].yflpar[j][1] >> theCurrentTemp.entx[k][i].yflpar[j][2]
387  >> theCurrentTemp.entx[k][i].yflpar[j][3] >> theCurrentTemp.entx[k][i].yflpar[j][4] >> theCurrentTemp.entx[k][i].yflpar[j][5];
388 
389  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 25, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
390  }
391 
392  for (j=0; j<4; ++j) {
393 
394  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];
395 
396  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 26, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
397  }
398 
399  for (j=0; j<4; ++j) {
400 
401  in_file >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >> theCurrentTemp.entx[k][i].xflpar[j][2]
402  >> theCurrentTemp.entx[k][i].xflpar[j][3] >> theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
403 
404  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 27, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
405  }
406 
407  for (j=0; j<4; ++j) {
408 
409  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];
410 
411  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 28, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
412  }
413 
414  for (j=0; j<4; ++j) {
415 
416  in_file >> theCurrentTemp.entx[k][i].yavgc2m[j] >> theCurrentTemp.entx[k][i].yrmsc2m[j] >> theCurrentTemp.entx[k][i].ygx0c2m[j] >> theCurrentTemp.entx[k][i].ygsigc2m[j];
417 
418  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 29, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
419  }
420 
421  for (j=0; j<4; ++j) {
422 
423  in_file >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >> theCurrentTemp.entx[k][i].xgx0c2m[j] >> theCurrentTemp.entx[k][i].xgsigc2m[j];
424 
425  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
426  }
427 
428  for (j=0; j<4; ++j) {
429 
430  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];
431 
432  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30a, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
433  }
434 
435  for (j=0; j<4; ++j) {
436 
437  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];
438 
439  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30b, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
440  }
441 
442  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
443  >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav >> theCurrentTemp.entx[k][i].qavg_spare >> theCurrentTemp.entx[k][i].spare[0];
444 
445  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 31, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
446 
447  in_file >> theCurrentTemp.entx[k][i].spare[1] >> theCurrentTemp.entx[k][i].spare[2] >> theCurrentTemp.entx[k][i].spare[3] >> theCurrentTemp.entx[k][i].qbfrac[0] >> theCurrentTemp.entx[k][i].qbfrac[1]
448  >> 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;
449  // 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];
450 
451  if(in_file.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 32, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
452 
453  }
454  }
455 
456 
457  in_file.close();
458 
459  // Add this template to the store
460 
461  thePixelTemp.push_back(theCurrentTemp);
462 
463  return true;
464 
465  } else {
466 
467  // If file didn't open, report this
468 
469  LOGERROR("SiPixelTemplate") << "Error opening File" << tempfile << ENDL;
470  return false;
471 
472  }
473 
474 } // TempInit
475 
476 
477 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
478 
479 //****************************************************************
483 //****************************************************************
484 bool SiPixelTemplate::pushfile(const SiPixelTemplateDBObject& dbobject)
485 {
486  // Add template stored in external dbobject to theTemplateStore
487 
488  // Local variables
489  int i, j, k, l;
490  float qavg_avg;
491  // const char *tempfile;
492  const int code_version={16};
493 
494  // We must create a new object because dbobject must be a const and our stream must not be
495  SiPixelTemplateDBObject db = dbobject;
496 
497  // Create a local template storage entry
498  SiPixelTemplateStore theCurrentTemp;
499 
500  // Fill the template storage for each template calibration stored in the db
501  for(int m=0; m<db.numOfTempl(); ++m)
502  {
503 
504  // Read-in a header string first and print it
505 
506  SiPixelTemplateDBObject::char2float temp;
507  for (i=0; i<20; ++i) {
508  temp.f = db.sVector()[db.index()];
509  theCurrentTemp.head.title[4*i] = temp.c[0];
510  theCurrentTemp.head.title[4*i+1] = temp.c[1];
511  theCurrentTemp.head.title[4*i+2] = temp.c[2];
512  theCurrentTemp.head.title[4*i+3] = temp.c[3];
513  db.incrementIndex(1);
514  }
515  theCurrentTemp.head.title[79] = '\0';
516  LOGINFO("SiPixelTemplate") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
517 
518  // next, the header information
519 
520  db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
521  >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
522  >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
523 
524  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file, no template load" << ENDL; return false;}
525 
526  LOGINFO("SiPixelTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
527  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
528  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
529  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
530  << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
531  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
532 
533  if(theCurrentTemp.head.templ_version < code_version) {LOGERROR("SiPixelTemplate") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
534 
535  // next, loop over all barrel y-angle entries
536 
537  for (i=0; i < theCurrentTemp.head.NTy; ++i) {
538 
539  db >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0]
540  >> theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2];
541 
542  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
543 
544  // Calculate the alpha, beta, and cot(beta) for this entry
545 
546  theCurrentTemp.enty[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
547 
548  theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0]/theCurrentTemp.enty[i].costrk[2];
549 
550  theCurrentTemp.enty[i].beta = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
551 
552  theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1]/theCurrentTemp.enty[i].costrk[2];
553 
554  db >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].symax >> theCurrentTemp.enty[i].dyone
555  >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].sxmax >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
556 
557  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
558 
559  db >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo
560  >> theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clsleny >> theCurrentTemp.enty[i].clslenx;
561  // >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav;
562 
563  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 3, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
564 
565  for (j=0; j<2; ++j) {
566 
567  db >> theCurrentTemp.enty[i].ypar[j][0] >> theCurrentTemp.enty[i].ypar[j][1]
568  >> theCurrentTemp.enty[i].ypar[j][2] >> theCurrentTemp.enty[i].ypar[j][3] >> theCurrentTemp.enty[i].ypar[j][4];
569 
570  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 4, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
571 
572  }
573 
574  for (j=0; j<9; ++j) {
575 
576  for (k=0; k<TYSIZE; ++k) {db >> theCurrentTemp.enty[i].ytemp[j][k];}
577 
578  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 5, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
579  }
580 
581  for (j=0; j<2; ++j) {
582 
583  db >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1]
584  >> theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
585 
586  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 6, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
587 
588  }
589 
590  qavg_avg = 0.;
591  for (j=0; j<9; ++j) {
592 
593  for (k=0; k<TXSIZE; ++k) {db >> theCurrentTemp.enty[i].xtemp[j][k]; qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];}
594 
595  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 7, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
596  }
597  theCurrentTemp.enty[i].qavg_avg = qavg_avg/9.;
598 
599  for (j=0; j<4; ++j) {
600 
601  db >> theCurrentTemp.enty[i].yavg[j] >> theCurrentTemp.enty[i].yrms[j] >> theCurrentTemp.enty[i].ygx0[j] >> theCurrentTemp.enty[i].ygsig[j];
602 
603  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 8, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
604  }
605 
606  for (j=0; j<4; ++j) {
607 
608  db >> theCurrentTemp.enty[i].yflpar[j][0] >> theCurrentTemp.enty[i].yflpar[j][1] >> theCurrentTemp.enty[i].yflpar[j][2]
609  >> theCurrentTemp.enty[i].yflpar[j][3] >> theCurrentTemp.enty[i].yflpar[j][4] >> theCurrentTemp.enty[i].yflpar[j][5];
610 
611  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 9, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
612  }
613 
614  for (j=0; j<4; ++j) {
615 
616  db >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >> theCurrentTemp.enty[i].xgsig[j];
617 
618  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 10, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
619  }
620 
621  for (j=0; j<4; ++j) {
622 
623  db >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >> theCurrentTemp.enty[i].xflpar[j][2]
624  >> theCurrentTemp.enty[i].xflpar[j][3] >> theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
625 
626  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 11, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
627  }
628 
629  for (j=0; j<4; ++j) {
630 
631  db >> theCurrentTemp.enty[i].chi2yavg[j] >> theCurrentTemp.enty[i].chi2ymin[j] >> theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j];
632 
633  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 12, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
634  }
635 
636  for (j=0; j<4; ++j) {
637 
638  db >> theCurrentTemp.enty[i].yavgc2m[j] >> theCurrentTemp.enty[i].yrmsc2m[j] >> theCurrentTemp.enty[i].ygx0c2m[j] >> theCurrentTemp.enty[i].ygsigc2m[j];
639 
640  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 13, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
641  }
642 
643  for (j=0; j<4; ++j) {
644 
645  db >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >> theCurrentTemp.enty[i].xgx0c2m[j] >> theCurrentTemp.enty[i].xgsigc2m[j];
646 
647  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
648  }
649 
650  for (j=0; j<4; ++j) {
651 
652  db >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >> theCurrentTemp.enty[i].ygx0gen[j] >> theCurrentTemp.enty[i].ygsiggen[j];
653 
654  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14a, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
655  }
656 
657  for (j=0; j<4; ++j) {
658 
659  db >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >> theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
660 
661  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 14b, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
662  }
663 
664 
665  db >> theCurrentTemp.enty[i].chi2yavgone >> theCurrentTemp.enty[i].chi2yminone >> theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2
666  >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav >> theCurrentTemp.enty[i].qavg_spare >> theCurrentTemp.enty[i].spare[0];
667 
668  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 15, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
669 
670  db >> theCurrentTemp.enty[i].spare[1] >> theCurrentTemp.enty[i].spare[2] >> theCurrentTemp.enty[i].spare[3] >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1]
671  >> theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracyone >> theCurrentTemp.enty[i].fracxone >> theCurrentTemp.enty[i].fracytwo >> theCurrentTemp.enty[i].fracxtwo;
672  // theCurrentTemp.enty[i].qbfrac[3] = 1. - theCurrentTemp.enty[i].qbfrac[0] - theCurrentTemp.enty[i].qbfrac[1] - theCurrentTemp.enty[i].qbfrac[2];
673 
674  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 16, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
675 
676  }
677 
678  // next, loop over all barrel x-angle entries
679 
680  for (k=0; k < theCurrentTemp.head.NTyx; ++k) {
681 
682  for (i=0; i < theCurrentTemp.head.NTxx; ++i) {
683 
684  db >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0]
685  >> theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2];
686 
687  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 17, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
688 
689  // Calculate the alpha, beta, and cot(beta) for this entry
690 
691  theCurrentTemp.entx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
692 
693  theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0]/theCurrentTemp.entx[k][i].costrk[2];
694 
695  theCurrentTemp.entx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
696 
697  theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1]/theCurrentTemp.entx[k][i].costrk[2];
698 
699  db >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >> theCurrentTemp.entx[k][i].symax >> theCurrentTemp.entx[k][i].dyone
700  >> theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
701 
702  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 18, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
703 
704  db >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >> theCurrentTemp.entx[k][i].dxtwo
705  >> theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clsleny >> theCurrentTemp.entx[k][i].clslenx;
706  // >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
707 
708  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 19, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
709 
710  for (j=0; j<2; ++j) {
711 
712  db >> theCurrentTemp.entx[k][i].ypar[j][0] >> theCurrentTemp.entx[k][i].ypar[j][1]
713  >> theCurrentTemp.entx[k][i].ypar[j][2] >> theCurrentTemp.entx[k][i].ypar[j][3] >> theCurrentTemp.entx[k][i].ypar[j][4];
714 
715  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 20, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
716  }
717 
718  for (j=0; j<9; ++j) {
719 
720  for (l=0; l<TYSIZE; ++l) {db >> theCurrentTemp.entx[k][i].ytemp[j][l];}
721 
722  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 21, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
723  }
724 
725  for (j=0; j<2; ++j) {
726 
727  db >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1]
728  >> theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >> theCurrentTemp.entx[k][i].xpar[j][4];
729 
730 
731  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 22, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
732  }
733 
734  qavg_avg = 0.;
735  for (j=0; j<9; ++j) {
736 
737  for (l=0; l<TXSIZE; ++l) {db >> theCurrentTemp.entx[k][i].xtemp[j][l]; qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];}
738 
739  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 23, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
740  }
741  theCurrentTemp.entx[k][i].qavg_avg = qavg_avg/9.;
742 
743  for (j=0; j<4; ++j) {
744 
745  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];
746 
747  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 24, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
748  }
749 
750  for (j=0; j<4; ++j) {
751 
752  db >> theCurrentTemp.entx[k][i].yflpar[j][0] >> theCurrentTemp.entx[k][i].yflpar[j][1] >> theCurrentTemp.entx[k][i].yflpar[j][2]
753  >> theCurrentTemp.entx[k][i].yflpar[j][3] >> theCurrentTemp.entx[k][i].yflpar[j][4] >> theCurrentTemp.entx[k][i].yflpar[j][5];
754 
755  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 25, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
756  }
757 
758  for (j=0; j<4; ++j) {
759 
760  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];
761 
762  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 26, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
763  }
764 
765  for (j=0; j<4; ++j) {
766 
767  db >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >> theCurrentTemp.entx[k][i].xflpar[j][2]
768  >> theCurrentTemp.entx[k][i].xflpar[j][3] >> theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
769 
770  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 27, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
771  }
772 
773  for (j=0; j<4; ++j) {
774 
775  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];
776 
777  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 28, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
778  }
779 
780  for (j=0; j<4; ++j) {
781 
782  db >> theCurrentTemp.entx[k][i].yavgc2m[j] >> theCurrentTemp.entx[k][i].yrmsc2m[j] >> theCurrentTemp.entx[k][i].ygx0c2m[j] >> theCurrentTemp.entx[k][i].ygsigc2m[j];
783 
784  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 29, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
785  }
786 
787  for (j=0; j<4; ++j) {
788 
789  db >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >> theCurrentTemp.entx[k][i].xgx0c2m[j] >> theCurrentTemp.entx[k][i].xgsigc2m[j];
790 
791  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
792  }
793 
794  for (j=0; j<4; ++j) {
795 
796  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];
797 
798  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30a, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
799  }
800 
801  for (j=0; j<4; ++j) {
802 
803  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];
804 
805  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 30b, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
806  }
807 
808 
809  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
810  >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav >> theCurrentTemp.entx[k][i].qavg_spare >> theCurrentTemp.entx[k][i].spare[0];
811 
812  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 31, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
813 
814  db >> theCurrentTemp.entx[k][i].spare[1] >> theCurrentTemp.entx[k][i].spare[2] >> theCurrentTemp.entx[k][i].spare[3] >> theCurrentTemp.entx[k][i].qbfrac[0] >> theCurrentTemp.entx[k][i].qbfrac[1]
815  >> 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;
816  // 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];
817 
818  if(db.fail()) {LOGERROR("SiPixelTemplate") << "Error reading file 32, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
819 
820  }
821  }
822 
823 
824  // Add this template to the store
825 
826  thePixelTemp.push_back(theCurrentTemp);
827 
828  }
829  return true;
830 
831 } // TempInit
832 
833 #endif
834 
835 
836 // ************************************************************************************************************
843 // ************************************************************************************************************
844 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBz) {
845  // Interpolate for a new set of track angles
846 
847  // Local variables
848  int i, j;
849  int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx;
850  float yratio, yxratio, xxratio, sxmax, qcorrect, qxtempcor, symax, chi2xavgone, chi2xminone, cotb, cotalpha0, cotbeta0;
851  bool flip_y;
852  // std::vector <float> xrms(4), xgsig(4), xrmsc2m(4), xgsigc2m(4);
853  float chi2xavg[4], chi2xmin[4];
854 
855 
856  // Check to see if interpolation is valid
857 
858  if(id != id_current || cotalpha != cota_current || cotbeta != cotb_current) {
859 
860  cota_current = cotalpha; cotb_current = cotbeta; success = true;
861 
862  if(id != id_current) {
863 
864  // Find the index corresponding to id
865 
866  index_id = -1;
867  for(i=0; i<(int)thePixelTemp.size(); ++i) {
868 
869  if(id == thePixelTemp[i].head.ID) {
870 
871  index_id = i;
872  id_current = id;
873 
874  // Copy the charge scaling factor to the private variable
875 
876  pqscale = thePixelTemp[index_id].head.qscale;
877 
878  // Copy the pseudopixel signal size to the private variable
879 
880  ps50 = thePixelTemp[index_id].head.s50;
881 
882  // Pixel sizes to the private variables
883 
884  pxsize = thePixelTemp[index_id].head.xsize;
885  pysize = thePixelTemp[index_id].head.ysize;
886  pzsize = thePixelTemp[index_id].head.zsize;
887 
888  break;
889  }
890  }
891  }
892 
893 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
894  if(index_id < 0 || index_id >= (int)thePixelTemp.size()) {
895  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::interpolate can't find needed template ID = " << id << std::endl;
896  }
897 #else
898  assert(index_id >= 0 && index_id < (int)thePixelTemp.size());
899 #endif
900 
901  // Interpolate the absolute value of cot(beta)
902 
903  abs_cotb = std::abs(cotbeta);
904 
905  // qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
906 
907  cotalpha0 = thePixelTemp[index_id].enty[0].cotalpha;
908  qcorrect= std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta*cotbeta+cotalpha0*cotalpha0));
909 
910 // for some cosmics, the ususal gymnastics are incorrect
911  if(thePixelTemp[index_id].head.Dtype == 0) {
912  cotb = abs_cotb;
913  flip_y = false;
914  if(cotbeta < 0.) {flip_y = true;}
915  } else {
916  if(locBz < 0.) {
917  cotb = cotbeta;
918  flip_y = false;
919  } else {
920  cotb = -cotbeta;
921  flip_y = true;
922  }
923  }
924 
925  Ny = thePixelTemp[index_id].head.NTy;
926  Nyx = thePixelTemp[index_id].head.NTyx;
927  Nxx = thePixelTemp[index_id].head.NTxx;
928 
929 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
930  if(Ny < 2 || Nyx < 1 || Nxx < 2) {
931  throw cms::Exception("DataCorrupt") << "template ID = " << id_current << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
932  }
933 #else
934  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
935 #endif
936  imaxx = Nyx - 1;
937  imidy = Nxx/2;
938 
939  // next, loop over all y-angle entries
940 
941  ilow = 0;
942  yratio = 0.;
943 
944  if(cotb >= thePixelTemp[index_id].enty[Ny-1].cotbeta) {
945 
946  ilow = Ny-2;
947  yratio = 1.;
948  success = false;
949 
950  } else {
951 
952  if(cotb >= thePixelTemp[index_id].enty[0].cotbeta) {
953 
954  for (i=0; i<Ny-1; ++i) {
955 
956  if( thePixelTemp[index_id].enty[i].cotbeta <= cotb && cotb < thePixelTemp[index_id].enty[i+1].cotbeta) {
957 
958  ilow = i;
959  yratio = (cotb - thePixelTemp[index_id].enty[i].cotbeta)/(thePixelTemp[index_id].enty[i+1].cotbeta - thePixelTemp[index_id].enty[i].cotbeta);
960  break;
961  }
962  }
963  } else { success = false; }
964  }
965 
966  ihigh=ilow + 1;
967 
968  // Interpolate/store all y-related quantities (flip displacements when flip_y)
969 
970  pyratio = yratio;
971  pqavg = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].qavg + yratio*thePixelTemp[index_id].enty[ihigh].qavg;
972  pqavg *= qcorrect;
973  symax = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].symax + yratio*thePixelTemp[index_id].enty[ihigh].symax;
974  psyparmax = symax;
975  sxmax = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].sxmax + yratio*thePixelTemp[index_id].enty[ihigh].sxmax;
976  pdyone = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].dyone + yratio*thePixelTemp[index_id].enty[ihigh].dyone;
977  if(flip_y) {pdyone = -pdyone;}
978  psyone = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].syone + yratio*thePixelTemp[index_id].enty[ihigh].syone;
979  pdytwo = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].dytwo + yratio*thePixelTemp[index_id].enty[ihigh].dytwo;
980  if(flip_y) {pdytwo = -pdytwo;}
981  psytwo = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].sytwo + yratio*thePixelTemp[index_id].enty[ihigh].sytwo;
982  pqmin = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].qmin + yratio*thePixelTemp[index_id].enty[ihigh].qmin;
983  pqmin *= qcorrect;
984  pqmin2 = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].qmin2 + yratio*thePixelTemp[index_id].enty[ihigh].qmin2;
985  pqmin2 *= qcorrect;
986  pmpvvav = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].mpvvav + yratio*thePixelTemp[index_id].enty[ihigh].mpvvav;
987  pmpvvav *= qcorrect;
988  psigmavav = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].sigmavav + yratio*thePixelTemp[index_id].enty[ihigh].sigmavav;
989  pkappavav = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].kappavav + yratio*thePixelTemp[index_id].enty[ihigh].kappavav;
990  pclsleny = fminf(thePixelTemp[index_id].enty[ilow].clsleny, thePixelTemp[index_id].enty[ihigh].clsleny);
991  pqavg_avg = (1.f - yratio)*thePixelTemp[index_id].enty[ilow].qavg_avg + yratio*thePixelTemp[index_id].enty[ihigh].qavg_avg;
992  pqavg_avg *= qcorrect;
993  for(i=0; i<2 ; ++i) {
994  for(j=0; j<5 ; ++j) {
995  // Charge loss switches sides when cot(beta) changes sign
996  if(flip_y) {
997  pyparl[1-i][j] = thePixelTemp[index_id].enty[ilow].ypar[i][j];
998  pyparh[1-i][j] = thePixelTemp[index_id].enty[ihigh].ypar[i][j];
999  } else {
1000  pyparl[i][j] = thePixelTemp[index_id].enty[ilow].ypar[i][j];
1001  pyparh[i][j] = thePixelTemp[index_id].enty[ihigh].ypar[i][j];
1002  }
1003  pxparly0[i][j] = thePixelTemp[index_id].enty[ilow].xpar[i][j];
1004  pxparhy0[i][j] = thePixelTemp[index_id].enty[ihigh].xpar[i][j];
1005  }
1006  }
1007  for(i=0; i<4; ++i) {
1008  pyavg[i]=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].yavg[i] + yratio*thePixelTemp[index_id].enty[ihigh].yavg[i];
1009  if(flip_y) {pyavg[i] = -pyavg[i];}
1010  pyrms[i]=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].yrms[i] + yratio*thePixelTemp[index_id].enty[ihigh].yrms[i];
1011  // pygx0[i]=(1. - yratio)*thePixelTemp[index_id].enty[ilow].ygx0[i] + yratio*thePixelTemp[index_id].enty[ihigh].ygx0[i];
1012  // if(flip_y) {pygx0[i] = -pygx0[i];}
1013  // pygsig[i]=(1. - yratio)*thePixelTemp[index_id].enty[ilow].ygsig[i] + yratio*thePixelTemp[index_id].enty[ihigh].ygsig[i];
1014  // xrms[i]=(1. - yratio)*thePixelTemp[index_id].enty[ilow].xrms[i] + yratio*thePixelTemp[index_id].enty[ihigh].xrms[i];
1015  // xgsig[i]=(1. - yratio)*thePixelTemp[index_id].enty[ilow].xgsig[i] + yratio*thePixelTemp[index_id].enty[ihigh].xgsig[i];
1016  pchi2yavg[i]=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].chi2yavg[i] + yratio*thePixelTemp[index_id].enty[ihigh].chi2yavg[i];
1017  pchi2ymin[i]=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].chi2ymin[i] + yratio*thePixelTemp[index_id].enty[ihigh].chi2ymin[i];
1018  chi2xavg[i]=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].chi2xavg[i] + yratio*thePixelTemp[index_id].enty[ihigh].chi2xavg[i];
1019  chi2xmin[i]=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].chi2xmin[i] + yratio*thePixelTemp[index_id].enty[ihigh].chi2xmin[i];
1020  pyavgc2m[i]=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].yavgc2m[i] + yratio*thePixelTemp[index_id].enty[ihigh].yavgc2m[i];
1021  if(flip_y) {pyavgc2m[i] = -pyavgc2m[i];}
1022  pyrmsc2m[i]=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].yrmsc2m[i] + yratio*thePixelTemp[index_id].enty[ihigh].yrmsc2m[i];
1023  // pygx0c2m[i]=(1. - yratio)*thePixelTemp[index_id].enty[ilow].ygx0c2m[i] + yratio*thePixelTemp[index_id].enty[ihigh].ygx0c2m[i];
1024  // if(flip_y) {pygx0c2m[i] = -pygx0c2m[i];}
1025  // pygsigc2m[i]=(1. - yratio)*thePixelTemp[index_id].enty[ilow].ygsigc2m[i] + yratio*thePixelTemp[index_id].enty[ihigh].ygsigc2m[i];
1026  // xrmsc2m[i]=(1. - yratio)*thePixelTemp[index_id].enty[ilow].xrmsc2m[i] + yratio*thePixelTemp[index_id].enty[ihigh].xrmsc2m[i];
1027  // xgsigc2m[i]=(1. - yratio)*thePixelTemp[index_id].enty[ilow].xgsigc2m[i] + yratio*thePixelTemp[index_id].enty[ihigh].xgsigc2m[i];
1028  for(j=0; j<6 ; ++j) {
1029  pyflparl[i][j] = thePixelTemp[index_id].enty[ilow].yflpar[i][j];
1030  pyflparh[i][j] = thePixelTemp[index_id].enty[ihigh].yflpar[i][j];
1031 
1032  // Since Q_fl is odd under cotbeta, it flips qutomatically, change only even terms
1033 
1034  if(flip_y && (j == 0 || j == 2 || j == 4)) {
1035  pyflparl[i][j] = - pyflparl[i][j];
1036  pyflparh[i][j] = - pyflparh[i][j];
1037  }
1038  }
1039  }
1040 
1042 
1043  pchi2yavgone=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].chi2yavgone + yratio*thePixelTemp[index_id].enty[ihigh].chi2yavgone;
1044  pchi2yminone=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].chi2yminone + yratio*thePixelTemp[index_id].enty[ihigh].chi2yminone;
1045  chi2xavgone=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].chi2xavgone + yratio*thePixelTemp[index_id].enty[ihigh].chi2xavgone;
1046  chi2xminone=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].chi2xminone + yratio*thePixelTemp[index_id].enty[ihigh].chi2xminone;
1047  // for(i=0; i<10; ++i) {
1048  // pyspare[i]=(1. - yratio)*thePixelTemp[index_id].enty[ilow].yspare[i] + yratio*thePixelTemp[index_id].enty[ihigh].yspare[i];
1049  // }
1050 
1051  // Interpolate and build the y-template
1052 
1053  for(i=0; i<9; ++i) {
1054  pytemp[i][0] = 0.f;
1055  pytemp[i][1] = 0.f;
1056  pytemp[i][BYM2] = 0.f;
1057  pytemp[i][BYM1] = 0.f;
1058  for(j=0; j<TYSIZE; ++j) {
1059 
1060  // Flip the basic y-template when the cotbeta is negative
1061 
1062  if(flip_y) {
1063  pytemp[8-i][BYM3-j]=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].ytemp[i][j] + yratio*thePixelTemp[index_id].enty[ihigh].ytemp[i][j];
1064  } else {
1065  pytemp[i][j+2]=(1.f - yratio)*thePixelTemp[index_id].enty[ilow].ytemp[i][j] + yratio*thePixelTemp[index_id].enty[ihigh].ytemp[i][j];
1066  }
1067  }
1068  }
1069 
1070  // next, loop over all x-angle entries, first, find relevant y-slices
1071 
1072  iylow = 0;
1073  yxratio = 0.f;
1074 
1075  if(abs_cotb >= thePixelTemp[index_id].entx[Nyx-1][0].cotbeta) {
1076 
1077  iylow = Nyx-2;
1078  yxratio = 1.f;
1079 
1080  } else if(abs_cotb >= thePixelTemp[index_id].entx[0][0].cotbeta) {
1081 
1082  for (i=0; i<Nyx-1; ++i) {
1083 
1084  if( thePixelTemp[index_id].entx[i][0].cotbeta <= abs_cotb && abs_cotb < thePixelTemp[index_id].entx[i+1][0].cotbeta) {
1085 
1086  iylow = i;
1087  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);
1088  break;
1089  }
1090  }
1091  }
1092 
1093  iyhigh=iylow + 1;
1094 
1095  ilow = 0;
1096  xxratio = 0.f;
1097 
1098  if(cotalpha >= thePixelTemp[index_id].entx[0][Nxx-1].cotalpha) {
1099 
1100  ilow = Nxx-2;
1101  xxratio = 1.f;
1102  success = false;
1103 
1104  } else {
1105 
1106  if(cotalpha >= thePixelTemp[index_id].entx[0][0].cotalpha) {
1107 
1108  for (i=0; i<Nxx-1; ++i) {
1109 
1110  if( thePixelTemp[index_id].entx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp[index_id].entx[0][i+1].cotalpha) {
1111 
1112  ilow = i;
1113  xxratio = (cotalpha - thePixelTemp[index_id].entx[0][i].cotalpha)/(thePixelTemp[index_id].entx[0][i+1].cotalpha - thePixelTemp[index_id].entx[0][i].cotalpha);
1114  break;
1115  }
1116  }
1117  } else { success = false; }
1118  }
1119 
1120  ihigh=ilow + 1;
1121 
1122  // Interpolate/store all x-related quantities
1123 
1124  pyxratio = yxratio;
1125  pxxratio = xxratio;
1126 
1127  // sxparmax defines the maximum charge for which the parameters xpar are defined (not rescaled by cotbeta)
1128 
1129  psxparmax = (1.f - xxratio)*thePixelTemp[index_id].entx[imaxx][ilow].sxmax + xxratio*thePixelTemp[index_id].entx[imaxx][ihigh].sxmax;
1130  psxmax = psxparmax;
1131  if(thePixelTemp[index_id].entx[imaxx][imidy].sxmax != 0.f) {psxmax=psxmax/thePixelTemp[index_id].entx[imaxx][imidy].sxmax*sxmax;}
1132  psymax = (1.f - xxratio)*thePixelTemp[index_id].entx[imaxx][ilow].symax + xxratio*thePixelTemp[index_id].entx[imaxx][ihigh].symax;
1133  if(thePixelTemp[index_id].entx[imaxx][imidy].symax != 0.) {psymax=psymax/thePixelTemp[index_id].entx[imaxx][imidy].symax*symax;}
1134  pdxone = (1.f - xxratio)*thePixelTemp[index_id].entx[0][ilow].dxone + xxratio*thePixelTemp[index_id].entx[0][ihigh].dxone;
1135  psxone = (1.f - xxratio)*thePixelTemp[index_id].entx[0][ilow].sxone + xxratio*thePixelTemp[index_id].entx[0][ihigh].sxone;
1136  pdxtwo = (1.f - xxratio)*thePixelTemp[index_id].entx[0][ilow].dxtwo + xxratio*thePixelTemp[index_id].entx[0][ihigh].dxtwo;
1137  psxtwo = (1.f - xxratio)*thePixelTemp[index_id].entx[0][ilow].sxtwo + xxratio*thePixelTemp[index_id].entx[0][ihigh].sxtwo;
1138  pclslenx = fminf(thePixelTemp[index_id].entx[0][ilow].clslenx, thePixelTemp[index_id].entx[0][ihigh].clslenx);
1139  for(i=0; i<2 ; ++i) {
1140  for(j=0; j<5 ; ++j) {
1141  pxpar0[i][j] = thePixelTemp[index_id].entx[imaxx][imidy].xpar[i][j];
1142  pxparl[i][j] = thePixelTemp[index_id].entx[imaxx][ilow].xpar[i][j];
1143  pxparh[i][j] = thePixelTemp[index_id].entx[imaxx][ihigh].xpar[i][j];
1144  }
1145  }
1146 
1147  // pixmax is the maximum allowed pixel charge (used for truncation)
1148 
1149  ppixmax=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp[index_id].entx[iylow][ilow].pixmax + xxratio*thePixelTemp[index_id].entx[iylow][ihigh].pixmax)
1150  +yxratio*((1.f - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].pixmax + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].pixmax);
1151 
1152  for(i=0; i<4; ++i) {
1153  pxavg[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp[index_id].entx[iylow][ilow].xavg[i] + xxratio*thePixelTemp[index_id].entx[iylow][ihigh].xavg[i])
1154  +yxratio*((1.f - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].xavg[i] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].xavg[i]);
1155 
1156  pxrms[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp[index_id].entx[iylow][ilow].xrms[i] + xxratio*thePixelTemp[index_id].entx[iylow][ihigh].xrms[i])
1157  +yxratio*((1.f - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].xrms[i] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].xrms[i]);
1158 
1159  // pxgx0[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entx[iylow][ilow].xgx0[i] + xxratio*thePixelTemp[index_id].entx[iylow][ihigh].xgx0[i])
1160  // +yxratio*((1. - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].xgx0[i] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].xgx0[i]);
1161 
1162  // pxgsig[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entx[iylow][ilow].xgsig[i] + xxratio*thePixelTemp[index_id].entx[iylow][ihigh].xgsig[i])
1163  // +yxratio*((1. - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].xgsig[i] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].xgsig[i]);
1164 
1165  pxavgc2m[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp[index_id].entx[iylow][ilow].xavgc2m[i] + xxratio*thePixelTemp[index_id].entx[iylow][ihigh].xavgc2m[i])
1166  +yxratio*((1.f - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].xavgc2m[i] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].xavgc2m[i]);
1167 
1168  pxrmsc2m[i]=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp[index_id].entx[iylow][ilow].xrmsc2m[i] + xxratio*thePixelTemp[index_id].entx[iylow][ihigh].xrmsc2m[i])
1169  +yxratio*((1.f - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].xrmsc2m[i] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].xrmsc2m[i]);
1170 
1171  // pxgx0c2m[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entx[iylow][ilow].xgx0c2m[i] + xxratio*thePixelTemp[index_id].entx[iylow][ihigh].xgx0c2m[i])
1172  // +yxratio*((1. - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].xgx0c2m[i] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].xgx0c2m[i]);
1173 
1174  // pxgsigc2m[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entx[iylow][ilow].xgsigc2m[i] + xxratio*thePixelTemp[index_id].entx[iylow][ihigh].xgsigc2m[i])
1175  // +yxratio*((1. - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].xgsigc2m[i] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].xgsigc2m[i]);
1176  //
1177  // Try new interpolation scheme
1178  //
1179  // pchi2xavg[i]=((1. - xxratio)*thePixelTemp[index_id].entx[imaxx][ilow].chi2xavg[i] + xxratio*thePixelTemp[index_id].entx[imaxx][ihigh].chi2xavg[i]);
1180  // if(thePixelTemp[index_id].entx[imaxx][imidy].chi2xavg[i] != 0.) {pchi2xavg[i]=pchi2xavg[i]/thePixelTemp[index_id].entx[imaxx][imidy].chi2xavg[i]*chi2xavg[i];}
1181  //
1182  // pchi2xmin[i]=((1. - xxratio)*thePixelTemp[index_id].entx[imaxx][ilow].chi2xmin[i] + xxratio*thePixelTemp[index_id].entx[imaxx][ihigh].chi2xmin[i]);
1183  // if(thePixelTemp[index_id].entx[imaxx][imidy].chi2xmin[i] != 0.) {pchi2xmin[i]=pchi2xmin[i]/thePixelTemp[index_id].entx[imaxx][imidy].chi2xmin[i]*chi2xmin[i];}
1184  //
1185  pchi2xavg[i]=((1.f - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].chi2xavg[i] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].chi2xavg[i]);
1186  if(thePixelTemp[index_id].entx[iyhigh][imidy].chi2xavg[i] != 0.f) {pchi2xavg[i]=pchi2xavg[i]/thePixelTemp[index_id].entx[iyhigh][imidy].chi2xavg[i]*chi2xavg[i];}
1187 
1188  pchi2xmin[i]=((1.f - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].chi2xmin[i] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].chi2xmin[i]);
1189  if(thePixelTemp[index_id].entx[iyhigh][imidy].chi2xmin[i] != 0.f) {pchi2xmin[i]=pchi2xmin[i]/thePixelTemp[index_id].entx[iyhigh][imidy].chi2xmin[i]*chi2xmin[i];}
1190 
1191  for(j=0; j<6 ; ++j) {
1192  pxflparll[i][j] = thePixelTemp[index_id].entx[iylow][ilow].xflpar[i][j];
1193  pxflparlh[i][j] = thePixelTemp[index_id].entx[iylow][ihigh].xflpar[i][j];
1194  pxflparhl[i][j] = thePixelTemp[index_id].entx[iyhigh][ilow].xflpar[i][j];
1195  pxflparhh[i][j] = thePixelTemp[index_id].entx[iyhigh][ihigh].xflpar[i][j];
1196  }
1197  }
1198 
1199  // Do the spares next
1200 
1201  pchi2xavgone=((1.f - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].chi2xavgone + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].chi2xavgone);
1202  if(thePixelTemp[index_id].entx[iyhigh][imidy].chi2xavgone != 0.f) {pchi2xavgone=pchi2xavgone/thePixelTemp[index_id].entx[iyhigh][imidy].chi2xavgone*chi2xavgone;}
1203 
1204  pchi2xminone=((1.f - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].chi2xminone + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].chi2xminone);
1205  if(thePixelTemp[index_id].entx[iyhigh][imidy].chi2xminone != 0.f) {pchi2xminone=pchi2xminone/thePixelTemp[index_id].entx[iyhigh][imidy].chi2xminone*chi2xminone;}
1206  // for(i=0; i<10; ++i) {
1207  // pxspare[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index_id].entx[iylow][ilow].xspare[i] + xxratio*thePixelTemp[index_id].entx[iylow][ihigh].xspare[i])
1208  // +yxratio*((1. - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].xspare[i] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].xspare[i]);
1209  // }
1210 
1211  // Interpolate and build the x-template
1212 
1213  // qxtempcor corrects the total charge to the actual track angles (not actually needed for the template fits, but useful for Guofan)
1214 
1215  cotbeta0 = thePixelTemp[index_id].entx[iyhigh][0].cotbeta;
1216  qxtempcor=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta0*cotbeta0+cotalpha*cotalpha));
1217 
1218  for(i=0; i<9; ++i) {
1219  pxtemp[i][0] = 0.f;
1220  pxtemp[i][1] = 0.f;
1221  pxtemp[i][BXM2] = 0.f;
1222  pxtemp[i][BXM1] = 0.f;
1223  for(j=0; j<TXSIZE; ++j) {
1224  // Take next largest x-slice for the x-template (it reduces bias in the forward direction after irradiation)
1225  // pxtemp[i][j+2]=(1. - xxratio)*thePixelTemp[index_id].entx[imaxx][ilow].xtemp[i][j] + xxratio*thePixelTemp[index_id].entx[imaxx][ihigh].xtemp[i][j];
1226  // pxtemp[i][j+2]=(1. - xxratio)*thePixelTemp[index_id].entx[iyhigh][ilow].xtemp[i][j] + xxratio*thePixelTemp[index_id].entx[iyhigh][ihigh].xtemp[i][j];
1227  pxtemp[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]);
1228  }
1229  }
1230 
1231  plorywidth = thePixelTemp[index_id].head.lorywidth;
1232  if(locBz > 0.f) {plorywidth = -plorywidth;}
1233  plorxwidth = thePixelTemp[index_id].head.lorxwidth;
1234 
1235  }
1236 
1237  return success;
1238 } // interpolate
1239 
1240 
1241 
1242 
1243 
1244 // ************************************************************************************************************
1249 // ************************************************************************************************************
1250 bool SiPixelTemplate::interpolate(int id, float cotalpha, float cotbeta)
1251 {
1252  // Interpolate for a new set of track angles
1253 
1254  // Local variables
1255  float locBz;
1256  locBz = -1.;
1257  if(cotbeta < 0.) {locBz = -locBz;}
1258  return SiPixelTemplate::interpolate(id, cotalpha, cotbeta, locBz);
1259 }
1260 
1261 
1262 
1263 
1264 // ************************************************************************************************************
1272 // ************************************************************************************************************
1273 void SiPixelTemplate::ysigma2(int fypix, int lypix, float sythr, float ysum[25], float ysig2[25]) {
1274  // Interpolate using quantities already stored in the private variables
1275 
1276  // Local variables
1277  int i;
1278  float sigi, sigi2, sigi3, sigi4, symax, qscale;
1279 
1280  // Make sure that input is OK
1281 
1282 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1283  if(fypix < 2 || fypix >= BYM2) {
1284  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with fypix = " << fypix << std::endl;
1285  }
1286 #else
1287  assert(fypix > 1 && fypix < BYM2);
1288 #endif
1289 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1290  if(lypix < fypix || lypix >= BYM2) {
1291  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with lypix/fypix = " << lypix << "/" << fypix << std::endl;
1292  }
1293 #else
1294  assert(lypix >= fypix && lypix < BYM2);
1295 #endif
1296 
1297  // Define the maximum signal to use in the parameterization
1298 
1299  symax = psymax;
1300  if(psymax > psyparmax) {symax = psyparmax;}
1301 
1302  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1303 
1304  for(i=fypix-2; i<=lypix+2; ++i) {
1305  if(i < fypix || i > lypix) {
1306 
1307  // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
1308 
1309  ysig2[i] = ps50*ps50;
1310  } else {
1311  if(ysum[i] < symax) {
1312  sigi = ysum[i];
1313  qscale = 1.;
1314  } else {
1315  sigi = symax;
1316  qscale = ysum[i]/symax;
1317  }
1318  sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
1319  if(i <= BHY) {
1320  ysig2[i] = (1.-pyratio)*
1321  (pyparl[0][0]+pyparl[0][1]*sigi+pyparl[0][2]*sigi2+pyparl[0][3]*sigi3+pyparl[0][4]*sigi4)
1322  + pyratio*
1323  (pyparh[0][0]+pyparh[0][1]*sigi+pyparh[0][2]*sigi2+pyparh[0][3]*sigi3+pyparh[0][4]*sigi4);
1324  } else {
1325  ysig2[i] = (1.-pyratio)*
1326  (pyparl[1][0]+pyparl[1][1]*sigi+pyparl[1][2]*sigi2+pyparl[1][3]*sigi3+pyparl[1][4]*sigi4)
1327  + pyratio*
1328  (pyparh[1][0]+pyparh[1][1]*sigi+pyparh[1][2]*sigi2+pyparh[1][3]*sigi3+pyparh[1][4]*sigi4);
1329  }
1330  ysig2[i] *=qscale;
1331  if(ysum[i] > sythr) {ysig2[i] = 1.e8;}
1332  if(ysig2[i] <= 0.) {LOGERROR("SiPixelTemplate") << "neg y-error-squared, id = " << id_current << ", index = " << index_id <<
1333  ", cot(alpha) = " << cota_current << ", cot(beta) = " << cotb_current << ", sigi = " << sigi << ENDL;}
1334  }
1335  }
1336 
1337  return;
1338 
1339 } // End ysigma2
1340 
1341 
1342 // ************************************************************************************************************
1348 // ************************************************************************************************************
1349 void SiPixelTemplate::ysigma2(float qpixel, int index, float& ysig2)
1350 
1351 {
1352  // Interpolate using quantities already stored in the private variables
1353 
1354  // Local variables
1355  float sigi, sigi2, sigi3, sigi4, symax, qscale, err2;
1356 
1357  // Make sure that input is OK
1358 
1359 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1360  if(index < 2 || index >= BYM2) {
1361  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ysigma2 called with index = " << index << std::endl;
1362  }
1363 #else
1364  assert(index > 1 && index < BYM2);
1365 #endif
1366 
1367  // Define the maximum signal to use in the parameterization
1368 
1369  symax = psymax;
1370  if(psymax > psyparmax) {symax = psyparmax;}
1371 
1372  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1373 
1374  if(qpixel < symax) {
1375  sigi = qpixel;
1376  qscale = 1.;
1377  } else {
1378  sigi = symax;
1379  qscale = qpixel/symax;
1380  }
1381  sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
1382  if(index <= BHY) {
1383  err2 = (1.-pyratio)*
1384  (pyparl[0][0]+pyparl[0][1]*sigi+pyparl[0][2]*sigi2+pyparl[0][3]*sigi3+pyparl[0][4]*sigi4)
1385  + pyratio*
1386  (pyparh[0][0]+pyparh[0][1]*sigi+pyparh[0][2]*sigi2+pyparh[0][3]*sigi3+pyparh[0][4]*sigi4);
1387  } else {
1388  err2 = (1.-pyratio)*
1389  (pyparl[1][0]+pyparl[1][1]*sigi+pyparl[1][2]*sigi2+pyparl[1][3]*sigi3+pyparl[1][4]*sigi4)
1390  + pyratio*
1391  (pyparh[1][0]+pyparh[1][1]*sigi+pyparh[1][2]*sigi2+pyparh[1][3]*sigi3+pyparh[1][4]*sigi4);
1392  }
1393  ysig2 =qscale*err2;
1394  if(ysig2 <= 0.) {LOGERROR("SiPixelTemplate") << "neg y-error-squared, id = " << id_current << ", index = " << index_id <<
1395  ", cot(alpha) = " << cota_current << ", cot(beta) = " << cotb_current << ", sigi = " << sigi << ENDL;}
1396 
1397  return;
1398 
1399 } // End ysigma2
1400 
1401 
1402 
1403 
1404 
1405 // ************************************************************************************************************
1413 // ************************************************************************************************************
1414  void SiPixelTemplate::xsigma2(int fxpix, int lxpix, float sxthr, float xsum[11], float xsig2[11]) {
1415  // Interpolate using quantities already stored in the private variables
1416 
1417  // Local variables
1418  int i;
1419  float sigi, sigi2, sigi3, sigi4, yint, sxmax, x0, qscale;
1420 
1421  // Make sure that input is OK
1422 
1423 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1424  if(fxpix < 2 || fxpix >= BXM2) {
1425  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xsigma2 called with fxpix = " << fxpix << std::endl;
1426  }
1427 #else
1428  assert(fxpix > 1 && fxpix < BXM2);
1429 #endif
1430 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1431  if(lxpix < fxpix || lxpix >= BXM2) {
1432  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xsigma2 called with lxpix/fxpix = " << lxpix << "/" << fxpix << std::endl;
1433  }
1434 #else
1435  assert(lxpix >= fxpix && lxpix < BXM2);
1436 #endif
1437 
1438 // Define the maximum signal to use in the parameterization
1439 
1440  sxmax = psxmax;
1441  if(psxmax > psxparmax) {sxmax = psxparmax;}
1442 
1443  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1444 
1445  for(i=fxpix-2; i<=lxpix+2; ++i) {
1446  if(i < fxpix || i > lxpix) {
1447 
1448  // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
1449 
1450  xsig2[i] = ps50*ps50;
1451  } else {
1452  if(xsum[i] < sxmax) {
1453  sigi = xsum[i];
1454  qscale = 1.;
1455  } else {
1456  sigi = sxmax;
1457  qscale = xsum[i]/sxmax;
1458  }
1459  sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
1460 
1461  // First, do the cotbeta interpolation
1462 
1463  if(i <= BHX) {
1464  yint = (1.-pyratio)*
1465  (pxparly0[0][0]+pxparly0[0][1]*sigi+pxparly0[0][2]*sigi2+pxparly0[0][3]*sigi3+pxparly0[0][4]*sigi4)
1466  + pyratio*
1467  (pxparhy0[0][0]+pxparhy0[0][1]*sigi+pxparhy0[0][2]*sigi2+pxparhy0[0][3]*sigi3+pxparhy0[0][4]*sigi4);
1468  } else {
1469  yint = (1.-pyratio)*
1470  (pxparly0[1][0]+pxparly0[1][1]*sigi+pxparly0[1][2]*sigi2+pxparly0[1][3]*sigi3+pxparly0[1][4]*sigi4)
1471  + pyratio*
1472  (pxparhy0[1][0]+pxparhy0[1][1]*sigi+pxparhy0[1][2]*sigi2+pxparhy0[1][3]*sigi3+pxparhy0[1][4]*sigi4);
1473  }
1474 
1475  // Next, do the cotalpha interpolation
1476 
1477  if(i <= BHX) {
1478  xsig2[i] = (1.-pxxratio)*
1479  (pxparl[0][0]+pxparl[0][1]*sigi+pxparl[0][2]*sigi2+pxparl[0][3]*sigi3+pxparl[0][4]*sigi4)
1480  + pxxratio*
1481  (pxparh[0][0]+pxparh[0][1]*sigi+pxparh[0][2]*sigi2+pxparh[0][3]*sigi3+pxparh[0][4]*sigi4);
1482  } else {
1483  xsig2[i] = (1.-pxxratio)*
1484  (pxparl[1][0]+pxparl[1][1]*sigi+pxparl[1][2]*sigi2+pxparl[1][3]*sigi3+pxparl[1][4]*sigi4)
1485  + pxxratio*
1486  (pxparh[1][0]+pxparh[1][1]*sigi+pxparh[1][2]*sigi2+pxparh[1][3]*sigi3+pxparh[1][4]*sigi4);
1487  }
1488 
1489  // Finally, get the mid-point value of the cotalpha function
1490 
1491  if(i <= BHX) {
1492  x0 = pxpar0[0][0]+pxpar0[0][1]*sigi+pxpar0[0][2]*sigi2+pxpar0[0][3]*sigi3+pxpar0[0][4]*sigi4;
1493  } else {
1494  x0 = pxpar0[1][0]+pxpar0[1][1]*sigi+pxpar0[1][2]*sigi2+pxpar0[1][3]*sigi3+pxpar0[1][4]*sigi4;
1495  }
1496 
1497  // Finally, rescale the yint value for cotalpha variation
1498 
1499  if(x0 != 0.) {xsig2[i] = xsig2[i]/x0 * yint;}
1500  xsig2[i] *=qscale;
1501  if(xsum[i] > sxthr) {xsig2[i] = 1.e8;}
1502  if(xsig2[i] <= 0.) {LOGERROR("SiPixelTemplate") << "neg x-error-squared, id = " << id_current << ", index = " << index_id <<
1503  ", cot(alpha) = " << cota_current << ", cot(beta) = " << cotb_current << ", sigi = " << sigi << ENDL;}
1504  }
1505  }
1506 
1507  return;
1508 
1509  } // End xsigma2
1510 
1511 
1512 
1513 
1514 
1515 
1516 // ************************************************************************************************************
1520 // ************************************************************************************************************
1521  float SiPixelTemplate::yflcorr(int binq, float qfly)
1522 
1523 {
1524  // Interpolate using quantities already stored in the private variables
1525 
1526  // Local variables
1527  float qfl, qfl2, qfl3, qfl4, qfl5, dy;
1528 
1529  // Make sure that input is OK
1530 
1531 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1532  if(binq < 0 || binq > 3) {
1533  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yflcorr called with binq = " << binq << std::endl;
1534  }
1535 #else
1536  assert(binq >= 0 && binq < 4);
1537 #endif
1538 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1539  if(fabs((double)qfly) > 1.) {
1540  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::yflcorr called with qfly = " << qfly << std::endl;
1541  }
1542 #else
1543  assert(fabs((double)qfly) <= 1.);
1544 #endif
1545 
1546 // Define the maximum signal to allow before de-weighting a pixel
1547 
1548  qfl = qfly;
1549 
1550  if(qfl < -0.9) {qfl = -0.9;}
1551  if(qfl > 0.9) {qfl = 0.9;}
1552 
1553 // Interpolate between the two polynomials
1554 
1555  qfl2 = qfl*qfl; qfl3 = qfl2*qfl; qfl4 = qfl3*qfl; qfl5 = qfl4*qfl;
1556  dy = (1.-pyratio)*(pyflparl[binq][0]+pyflparl[binq][1]*qfl+pyflparl[binq][2]*qfl2+pyflparl[binq][3]*qfl3+pyflparl[binq][4]*qfl4+pyflparl[binq][5]*qfl5)
1557  + pyratio*(pyflparh[binq][0]+pyflparh[binq][1]*qfl+pyflparh[binq][2]*qfl2+pyflparh[binq][3]*qfl3+pyflparh[binq][4]*qfl4+pyflparh[binq][5]*qfl5);
1558 
1559  return dy;
1560 
1561 } // End yflcorr
1562 
1563 
1564 
1565 
1566 
1567 
1568 // ************************************************************************************************************
1572 // ************************************************************************************************************
1573  float SiPixelTemplate::xflcorr(int binq, float qflx)
1574 
1575 {
1576  // Interpolate using quantities already stored in the private variables
1577 
1578  // Local variables
1579  float qfl, qfl2, qfl3, qfl4, qfl5, dx;
1580 
1581  // Make sure that input is OK
1582 
1583 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1584  if(binq < 0 || binq > 3) {
1585  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xflcorr called with binq = " << binq << std::endl;
1586  }
1587 #else
1588  assert(binq >= 0 && binq < 4);
1589 #endif
1590 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1591  if(fabs((double)qflx) > 1.) {
1592  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xflcorr called with qflx = " << qflx << std::endl;
1593  }
1594 #else
1595  assert(fabs((double)qflx) <= 1.);
1596 #endif
1597 
1598 // Define the maximum signal to allow before de-weighting a pixel
1599 
1600  qfl = qflx;
1601 
1602  if(qfl < -0.9) {qfl = -0.9;}
1603  if(qfl > 0.9) {qfl = 0.9;}
1604 
1605 // Interpolate between the two polynomials
1606 
1607  qfl2 = qfl*qfl; qfl3 = qfl2*qfl; qfl4 = qfl3*qfl; qfl5 = qfl4*qfl;
1608  dx = (1. - pyxratio)*((1.-pxxratio)*(pxflparll[binq][0]+pxflparll[binq][1]*qfl+pxflparll[binq][2]*qfl2+pxflparll[binq][3]*qfl3+pxflparll[binq][4]*qfl4+pxflparll[binq][5]*qfl5)
1609  + pxxratio*(pxflparlh[binq][0]+pxflparlh[binq][1]*qfl+pxflparlh[binq][2]*qfl2+pxflparlh[binq][3]*qfl3+pxflparlh[binq][4]*qfl4+pxflparlh[binq][5]*qfl5))
1610  + pyxratio*((1.-pxxratio)*(pxflparhl[binq][0]+pxflparhl[binq][1]*qfl+pxflparhl[binq][2]*qfl2+pxflparhl[binq][3]*qfl3+pxflparhl[binq][4]*qfl4+pxflparhl[binq][5]*qfl5)
1611  + pxxratio*(pxflparhh[binq][0]+pxflparhh[binq][1]*qfl+pxflparhh[binq][2]*qfl2+pxflparhh[binq][3]*qfl3+pxflparhh[binq][4]*qfl4+pxflparhh[binq][5]*qfl5));
1612 
1613  return dx;
1614 
1615 } // End xflcorr
1616 
1617 
1618 
1619 // ************************************************************************************************************
1624 // ************************************************************************************************************
1625  void SiPixelTemplate::ytemp(int fybin, int lybin, float ytemplate[41][BYSIZE])
1626 
1627 {
1628  // Retrieve already interpolated quantities
1629 
1630  // Local variables
1631  int i, j;
1632 
1633  // Verify that input parameters are in valid range
1634 
1635 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1636  if(fybin < 0 || fybin > 40) {
1637  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp called with fybin = " << fybin << std::endl;
1638  }
1639 #else
1640  assert(fybin >= 0 && fybin < 41);
1641 #endif
1642 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1643  if(lybin < 0 || lybin > 40) {
1644  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp called with lybin = " << lybin << std::endl;
1645  }
1646 #else
1647  assert(lybin >= 0 && lybin < 41);
1648 #endif
1649 
1650 // Build the y-template, the central 25 bins are here in all cases
1651 
1652  for(i=0; i<9; ++i) {
1653  for(j=0; j<BYSIZE; ++j) {
1654  ytemplate[i+16][j]=pytemp[i][j];
1655  }
1656  }
1657  for(i=0; i<8; ++i) {
1658  ytemplate[i+8][BYM1] = 0.;
1659  for(j=0; j<BYM1; ++j) {
1660  ytemplate[i+8][j]=pytemp[i][j+1];
1661  }
1662  }
1663  for(i=1; i<9; ++i) {
1664  ytemplate[i+24][0] = 0.;
1665  for(j=0; j<BYM1; ++j) {
1666  ytemplate[i+24][j+1]=pytemp[i][j];
1667  }
1668  }
1669 
1670 // Add more bins if needed
1671 
1672  if(fybin < 8) {
1673  for(i=0; i<8; ++i) {
1674  ytemplate[i][BYM2] = 0.;
1675  ytemplate[i][BYM1] = 0.;
1676  for(j=0; j<BYM2; ++j) {
1677  ytemplate[i][j]=pytemp[i][j+2];
1678  }
1679  }
1680  }
1681  if(lybin > 32) {
1682  for(i=1; i<9; ++i) {
1683  ytemplate[i+32][0] = 0.;
1684  ytemplate[i+32][1] = 0.;
1685  for(j=0; j<BYM2; ++j) {
1686  ytemplate[i+32][j+2]=pytemp[i][j];
1687  }
1688  }
1689  }
1690 
1691  return;
1692 
1693 } // End ytemp
1694 
1695 
1696 
1697 // ************************************************************************************************************
1702 // ************************************************************************************************************
1703  void SiPixelTemplate::xtemp(int fxbin, int lxbin, float xtemplate[41][BXSIZE])
1704 
1705 {
1706  // Retrieve already interpolated quantities
1707 
1708  // Local variables
1709  int i, j;
1710 
1711  // Verify that input parameters are in valid range
1712 
1713 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1714  if(fxbin < 0 || fxbin > 40) {
1715  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp called with fxbin = " << fxbin << std::endl;
1716  }
1717 #else
1718  assert(fxbin >= 0 && fxbin < 41);
1719 #endif
1720 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1721  if(lxbin < 0 || lxbin > 40) {
1722  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp called with lxbin = " << lxbin << std::endl;
1723  }
1724 #else
1725  assert(lxbin >= 0 && lxbin < 41);
1726 #endif
1727 
1728 // Build the x-template, the central 25 bins are here in all cases
1729 
1730  for(i=0; i<9; ++i) {
1731  for(j=0; j<BXSIZE; ++j) {
1732  xtemplate[i+16][j]=pxtemp[i][j];
1733  }
1734  }
1735  for(i=0; i<8; ++i) {
1736  xtemplate[i+8][BXM1] = 0.;
1737  for(j=0; j<BXM1; ++j) {
1738  xtemplate[i+8][j]=pxtemp[i][j+1];
1739  }
1740  }
1741  for(i=1; i<9; ++i) {
1742  xtemplate[i+24][0] = 0.;
1743  for(j=0; j<BXM1; ++j) {
1744  xtemplate[i+24][j+1]=pxtemp[i][j];
1745  }
1746  }
1747 
1748 // Add more bins if needed
1749 
1750  if(fxbin < 8) {
1751  for(i=0; i<8; ++i) {
1752  xtemplate[i][BXM2] = 0.;
1753  xtemplate[i][BXM1] = 0.;
1754  for(j=0; j<BXM2; ++j) {
1755  xtemplate[i][j]=pxtemp[i][j+2];
1756  }
1757  }
1758  }
1759  if(lxbin > 32) {
1760  for(i=1; i<9; ++i) {
1761  xtemplate[i+32][0] = 0.;
1762  xtemplate[i+32][1] = 0.;
1763  for(j=0; j<BXM2; ++j) {
1764  xtemplate[i+32][j+2]=pxtemp[i][j];
1765  }
1766  }
1767  }
1768 
1769  return;
1770 
1771 } // End xtemp
1772 
1773 
1774 // ************************************************************************************************************
1778 // ************************************************************************************************************
1779  void SiPixelTemplate::ytemp3d(int nypix, array_3d& ytemplate)
1780 
1781 {
1782  typedef boost::multi_array<float, 2> array_2d;
1783 
1784  // Retrieve already interpolated quantities
1785 
1786  // Local variables
1787  int i, j, k;
1788  int ioff0, ioffp, ioffm;
1789 
1790  // Verify that input parameters are in valid range
1791 
1792 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1793  if(nypix < 1 || nypix >= BYM3) {
1794  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::ytemp3d called with nypix = " << nypix << std::endl;
1795  }
1796 #else
1797  assert(nypix > 0 && nypix < BYM3);
1798 #endif
1799 
1800 // Calculate the size of the shift in pixels needed to span the entire cluster
1801 
1802  float diff = fabsf(nypix - pclsleny)/2. + 1.;
1803  int nshift = (int)diff;
1804  if((diff - nshift) > 0.5) {++nshift;}
1805 
1806 // Calculate the number of bins needed to specify each hit range
1807 
1808  int nbins = 9 + 16*nshift;
1809 
1810 // Create a 2-d working template with the correct size
1811 
1812  array_2d temp2d(boost::extents[nbins][BYSIZE]);
1813 
1814 // The 9 central bins are copied from the interpolated private store
1815 
1816  ioff0 = 8*nshift;
1817 
1818  for(i=0; i<9; ++i) {
1819  for(j=0; j<BYSIZE; ++j) {
1820  temp2d[i+ioff0][j]=pytemp[i][j];
1821  }
1822  }
1823 
1824 // Add the +- shifted templates
1825 
1826  for(k=1; k<=nshift; ++k) {
1827  ioffm=ioff0-k*8;
1828  for(i=0; i<8; ++i) {
1829  for(j=0; j<k; ++j) {
1830  temp2d[i+ioffm][BYM1-j] = 0.;
1831  }
1832  for(j=0; j<BYSIZE-k; ++j) {
1833  temp2d[i+ioffm][j]=pytemp[i][j+k];
1834  }
1835  }
1836  ioffp=ioff0+k*8;
1837  for(i=1; i<9; ++i) {
1838  for(j=0; j<k; ++j) {
1839  temp2d[i+ioffp][j] = 0.;
1840  }
1841  for(j=0; j<BYSIZE-k; ++j) {
1842  temp2d[i+ioffp][j+k]=pytemp[i][j];
1843  }
1844  }
1845  }
1846 
1847 // Resize the 3d template container
1848 
1849  ytemplate.resize(boost::extents[nbins][nbins][BYSIZE]);
1850 
1851 // Sum two 2-d templates to make the 3-d template
1852 
1853  for(i=0; i<nbins; ++i) {
1854  for(j=0; j<=i; ++j) {
1855  for(k=0; k<BYSIZE; ++k) {
1856  ytemplate[i][j][k]=temp2d[i][k]+temp2d[j][k];
1857  }
1858  }
1859  }
1860 
1861  return;
1862 
1863 } // End ytemp3d
1864 
1865 
1866 // ************************************************************************************************************
1870 // ************************************************************************************************************
1871  void SiPixelTemplate::xtemp3d(int nxpix, array_3d& xtemplate)
1872 
1873 {
1874  typedef boost::multi_array<float, 2> array_2d;
1875 
1876  // Retrieve already interpolated quantities
1877 
1878  // Local variables
1879  int i, j, k;
1880  int ioff0, ioffp, ioffm;
1881 
1882  // Verify that input parameters are in valid range
1883 
1884 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1885  if(nxpix < 1 || nxpix >= BXM3) {
1886  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp3d called with nxpix = " << nxpix << std::endl;
1887  }
1888 #else
1889  assert(nxpix > 0 && nxpix < BXM3);
1890 #endif
1891 
1892 // Calculate the size of the shift in pixels needed to span the entire cluster
1893 
1894  float diff = fabsf(nxpix - pclslenx)/2. + 1.;
1895  int nshift = (int)diff;
1896  if((diff - nshift) > 0.5) {++nshift;}
1897 
1898 // Calculate the number of bins needed to specify each hit range
1899 
1900  int nbins = 9 + 16*nshift;
1901 
1902 // Create a 2-d working template with the correct size
1903 
1904  array_2d temp2d(boost::extents[nbins][BXSIZE]);
1905 
1906 // The 9 central bins are copied from the interpolated private store
1907 
1908  ioff0 = 8*nshift;
1909 
1910  for(i=0; i<9; ++i) {
1911  for(j=0; j<BXSIZE; ++j) {
1912  temp2d[i+ioff0][j]=pxtemp[i][j];
1913  }
1914  }
1915 
1916 // Add the +- shifted templates
1917 
1918  for(k=1; k<=nshift; ++k) {
1919  ioffm=ioff0-k*8;
1920  for(i=0; i<8; ++i) {
1921  for(j=0; j<k; ++j) {
1922  temp2d[i+ioffm][BXM1-j] = 0.;
1923  }
1924  for(j=0; j<BXSIZE-k; ++j) {
1925  temp2d[i+ioffm][j]=pxtemp[i][j+k];
1926  }
1927  }
1928  ioffp=ioff0+k*8;
1929  for(i=1; i<9; ++i) {
1930  for(j=0; j<k; ++j) {
1931  temp2d[i+ioffp][j] = 0.;
1932  }
1933  for(j=0; j<BXSIZE-k; ++j) {
1934  temp2d[i+ioffp][j+k]=pxtemp[i][j];
1935  }
1936  }
1937  }
1938 
1939 // Resize the 3d template container
1940 
1941  xtemplate.resize(boost::extents[nbins][nbins][BXSIZE]);
1942 
1943 // Sum two 2-d templates to make the 3-d template
1944 
1945  for(i=0; i<nbins; ++i) {
1946  for(j=0; j<=i; ++j) {
1947  for(k=0; k<BXSIZE; ++k) {
1948  xtemplate[i][j][k]=temp2d[i][k]+temp2d[j][k];
1949  }
1950  }
1951  }
1952 
1953  return;
1954 
1955 } // End xtemp3d
1956 
1957 
1958 
1959 // ************************************************************************************************************
1983 // ************************************************************************************************************
1984 int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float& pixmx, float& sigmay, float& deltay, float& sigmax, float& deltax,
1985  float& sy1, float& dy1, float& sy2, float& dy2, float& sx1, float& dx1, float& sx2, float& dx2, float& lorywidth, float& lorxwidth)
1986 
1987 {
1988  // Interpolate for a new set of track angles
1989 
1990  // Local variables
1991  int i, binq;
1992  int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx, index;
1993  float yratio, yxratio, xxratio;
1994  float acotb, qscale, qavg, qmin, qmin2, fq, qtotal, qcorrect, cotb, cotalpha0;
1995  float yavggen[4], yrmsgen[4], xavggen[4], xrmsgen[4];
1996  bool flip_y;
1997 
1998 
1999 // Find the index corresponding to id
2000 
2001  index = -1;
2002  for(i=0; i<(int)thePixelTemp.size(); ++i) {
2003 
2004  if(id == thePixelTemp[i].head.ID) {
2005 
2006  index = i;
2007  break;
2008  }
2009  }
2010 
2011 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2012  if(index < 0 || index >= (int)thePixelTemp.size()) {
2013  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::qbin can't find needed template ID = " << id << std::endl;
2014  }
2015 #else
2016  assert(index >= 0 && index < (int)thePixelTemp.size());
2017 #endif
2018 
2019 //
2020 
2021 // Interpolate the absolute value of cot(beta)
2022 
2023  acotb = fabs((double)cotbeta);
2024 
2025 // qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
2026 
2027  // qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
2028 
2029  cotalpha0 = thePixelTemp[index].enty[0].cotalpha;
2030  qcorrect=(float)sqrt((double)((1.+cotbeta*cotbeta+cotalpha*cotalpha)/(1.+cotbeta*cotbeta+cotalpha0*cotalpha0)));
2031 
2032  // for some cosmics, the ususal gymnastics are incorrect
2033 
2034  if(thePixelTemp[index].head.Dtype == 0) {
2035  cotb = acotb;
2036  flip_y = false;
2037  if(cotbeta < 0.) {flip_y = true;}
2038  } else {
2039  if(locBz < 0.) {
2040  cotb = cotbeta;
2041  flip_y = false;
2042  } else {
2043  cotb = -cotbeta;
2044  flip_y = true;
2045  }
2046  }
2047 
2048  // Copy the charge scaling factor to the private variable
2049 
2050  qscale = thePixelTemp[index].head.qscale;
2051 
2052  Ny = thePixelTemp[index].head.NTy;
2053  Nyx = thePixelTemp[index].head.NTyx;
2054  Nxx = thePixelTemp[index].head.NTxx;
2055 
2056 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2057  if(Ny < 2 || Nyx < 1 || Nxx < 2) {
2058  throw cms::Exception("DataCorrupt") << "template ID = " << id_current << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
2059  }
2060 #else
2061  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
2062 #endif
2063  imaxx = Nyx - 1;
2064  imidy = Nxx/2;
2065 
2066 // next, loop over all y-angle entries
2067 
2068  ilow = 0;
2069  yratio = 0.;
2070 
2071  if(cotb >= thePixelTemp[index].enty[Ny-1].cotbeta) {
2072 
2073  ilow = Ny-2;
2074  yratio = 1.;
2075 
2076  } else {
2077 
2078  if(cotb >= thePixelTemp[index].enty[0].cotbeta) {
2079 
2080  for (i=0; i<Ny-1; ++i) {
2081 
2082  if( thePixelTemp[index].enty[i].cotbeta <= cotb && cotb < thePixelTemp[index].enty[i+1].cotbeta) {
2083 
2084  ilow = i;
2085  yratio = (cotb - thePixelTemp[index].enty[i].cotbeta)/(thePixelTemp[index].enty[i+1].cotbeta - thePixelTemp[index].enty[i].cotbeta);
2086  break;
2087  }
2088  }
2089  }
2090  }
2091 
2092  ihigh=ilow + 1;
2093 
2094 // Interpolate/store all y-related quantities (flip displacements when flip_y)
2095 
2096  qavg = (1. - yratio)*thePixelTemp[index].enty[ilow].qavg + yratio*thePixelTemp[index].enty[ihigh].qavg;
2097  qavg *= qcorrect;
2098  dy1 = (1. - yratio)*thePixelTemp[index].enty[ilow].dyone + yratio*thePixelTemp[index].enty[ihigh].dyone;
2099  if(flip_y) {dy1 = -dy1;}
2100  sy1 = (1. - yratio)*thePixelTemp[index].enty[ilow].syone + yratio*thePixelTemp[index].enty[ihigh].syone;
2101  dy2 = (1. - yratio)*thePixelTemp[index].enty[ilow].dytwo + yratio*thePixelTemp[index].enty[ihigh].dytwo;
2102  if(flip_y) {dy2 = -dy2;}
2103  sy2 = (1. - yratio)*thePixelTemp[index].enty[ilow].sytwo + yratio*thePixelTemp[index].enty[ihigh].sytwo;
2104  qmin = (1. - yratio)*thePixelTemp[index].enty[ilow].qmin + yratio*thePixelTemp[index].enty[ihigh].qmin;
2105  qmin *= qcorrect;
2106  qmin2 = (1. - yratio)*thePixelTemp[index].enty[ilow].qmin2 + yratio*thePixelTemp[index].enty[ihigh].qmin2;
2107  qmin2 *= qcorrect;
2108  for(i=0; i<4; ++i) {
2109  yavggen[i]=(1. - yratio)*thePixelTemp[index].enty[ilow].yavggen[i] + yratio*thePixelTemp[index].enty[ihigh].yavggen[i];
2110  if(flip_y) {yavggen[i] = -yavggen[i];}
2111  yrmsgen[i]=(1. - yratio)*thePixelTemp[index].enty[ilow].yrmsgen[i] + yratio*thePixelTemp[index].enty[ihigh].yrmsgen[i];
2112  }
2113 
2114 
2115 // next, loop over all x-angle entries, first, find relevant y-slices
2116 
2117  iylow = 0;
2118  yxratio = 0.;
2119 
2120  if(acotb >= thePixelTemp[index].entx[Nyx-1][0].cotbeta) {
2121 
2122  iylow = Nyx-2;
2123  yxratio = 1.;
2124 
2125  } else if(acotb >= thePixelTemp[index].entx[0][0].cotbeta) {
2126 
2127  for (i=0; i<Nyx-1; ++i) {
2128 
2129  if( thePixelTemp[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp[index].entx[i+1][0].cotbeta) {
2130 
2131  iylow = i;
2132  yxratio = (acotb - thePixelTemp[index].entx[i][0].cotbeta)/(thePixelTemp[index].entx[i+1][0].cotbeta - thePixelTemp[index].entx[i][0].cotbeta);
2133  break;
2134  }
2135  }
2136  }
2137 
2138  iyhigh=iylow + 1;
2139 
2140  ilow = 0;
2141  xxratio = 0.;
2142 
2143  if(cotalpha >= thePixelTemp[index].entx[0][Nxx-1].cotalpha) {
2144 
2145  ilow = Nxx-2;
2146  xxratio = 1.;
2147 
2148  } else {
2149 
2150  if(cotalpha >= thePixelTemp[index].entx[0][0].cotalpha) {
2151 
2152  for (i=0; i<Nxx-1; ++i) {
2153 
2154  if( thePixelTemp[index].entx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp[index].entx[0][i+1].cotalpha) {
2155 
2156  ilow = i;
2157  xxratio = (cotalpha - thePixelTemp[index].entx[0][i].cotalpha)/(thePixelTemp[index].entx[0][i+1].cotalpha - thePixelTemp[index].entx[0][i].cotalpha);
2158  break;
2159  }
2160  }
2161  }
2162  }
2163 
2164  ihigh=ilow + 1;
2165 
2166  dx1 = (1. - xxratio)*thePixelTemp[index].entx[0][ilow].dxone + xxratio*thePixelTemp[index].entx[0][ihigh].dxone;
2167  sx1 = (1. - xxratio)*thePixelTemp[index].entx[0][ilow].sxone + xxratio*thePixelTemp[index].entx[0][ihigh].sxone;
2168  dx2 = (1. - xxratio)*thePixelTemp[index].entx[0][ilow].dxtwo + xxratio*thePixelTemp[index].entx[0][ihigh].dxtwo;
2169  sx2 = (1. - xxratio)*thePixelTemp[index].entx[0][ilow].sxtwo + xxratio*thePixelTemp[index].entx[0][ihigh].sxtwo;
2170 
2171 // pixmax is the maximum allowed pixel charge (used for truncation)
2172 
2173  pixmx=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index].entx[iylow][ilow].pixmax + xxratio*thePixelTemp[index].entx[iylow][ihigh].pixmax)
2174  +yxratio*((1. - xxratio)*thePixelTemp[index].entx[iyhigh][ilow].pixmax + xxratio*thePixelTemp[index].entx[iyhigh][ihigh].pixmax);
2175 
2176  for(i=0; i<4; ++i) {
2177 
2178  xavggen[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index].entx[iylow][ilow].xavggen[i] + xxratio*thePixelTemp[index].entx[iylow][ihigh].xavggen[i])
2179  +yxratio*((1. - xxratio)*thePixelTemp[index].entx[iyhigh][ilow].xavggen[i] + xxratio*thePixelTemp[index].entx[iyhigh][ihigh].xavggen[i]);
2180 
2181  xrmsgen[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index].entx[iylow][ilow].xrmsgen[i] + xxratio*thePixelTemp[index].entx[iylow][ihigh].xrmsgen[i])
2182  +yxratio*((1. - xxratio)*thePixelTemp[index].entx[iyhigh][ilow].xrmsgen[i] + xxratio*thePixelTemp[index].entx[iyhigh][ihigh].xrmsgen[i]);
2183  }
2184 
2185  lorywidth = thePixelTemp[index].head.lorywidth;
2186  if(locBz > 0.) {lorywidth = -lorywidth;}
2187  lorxwidth = thePixelTemp[index].head.lorxwidth;
2188 
2189 
2190 
2191 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2192  if(qavg <= 0. || qmin <= 0.) {
2193  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::qbin, qavg or qmin <= 0,"
2194  << " Probably someone called the generic pixel reconstruction with an illegal trajectory state" << std::endl;
2195  }
2196 #else
2197  assert(qavg > 0. && qmin > 0.);
2198 #endif
2199 
2200 // Scale the input charge to account for differences between pixelav and CMSSW simulation or data
2201 
2202  qtotal = qscale*qclus;
2203 
2204 // uncertainty and final corrections depend upon total charge bin
2205 
2206  fq = qtotal/qavg;
2207  if(fq > 1.5) {
2208  binq=0;
2209  } else {
2210  if(fq > 1.0) {
2211  binq=1;
2212  } else {
2213  if(fq > 0.85) {
2214  binq=2;
2215  } else {
2216  binq=3;
2217  }
2218  }
2219  }
2220 
2221 // Take the errors and bias from the correct charge bin
2222 
2223  sigmay = yrmsgen[binq]; deltay = yavggen[binq];
2224 
2225  sigmax = xrmsgen[binq]; deltax = xavggen[binq];
2226 
2227 // If the charge is too small (then flag it)
2228 
2229  if(qtotal < 0.95*qmin) {binq = 5;} else {if(qtotal < 0.95*qmin2) {binq = 4;}}
2230 
2231  return binq;
2232 
2233 } // qbin
2234 
2235 // ************************************************************************************************************
2257 // ************************************************************************************************************
2258 int SiPixelTemplate::qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float& pixmx, float& sigmay, float& deltay, float& sigmax, float& deltax,
2259  float& sy1, float& dy1, float& sy2, float& dy2, float& sx1, float& dx1, float& sx2, float& dx2)
2260 
2261 {
2262  float lorywidth, lorxwidth;
2263  return SiPixelTemplate::qbin(id, cotalpha, cotbeta, locBz, qclus, pixmx, sigmay, deltay, sigmax, deltax,
2264  sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, lorywidth, lorxwidth);
2265 
2266 } // qbin
2267 
2268 // ************************************************************************************************************
2274 // ************************************************************************************************************
2275 int SiPixelTemplate::qbin(int id, float cotbeta, float qclus)
2276 {
2277 // Interpolate for a new set of track angles
2278 
2279 // Local variables
2280  float pixmx, sigmay, deltay, sigmax, deltax, sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, locBz, lorywidth, lorxwidth;
2281  const float cotalpha = 0.;
2282  locBz = -1.;
2283  if(cotbeta < 0.) {locBz = -locBz;}
2284  return SiPixelTemplate::qbin(id, cotalpha, cotbeta, locBz, qclus, pixmx, sigmay, deltay, sigmax, deltax,
2285  sy1, dy1, sy2, dy2, sx1, dx1, sx2, dx2, lorywidth, lorxwidth);
2286 
2287 } // qbin
2288 
2289 
2290 
2291 // ************************************************************************************************************
2303 // ************************************************************************************************************
2304 void SiPixelTemplate::temperrors(int id, float cotalpha, float cotbeta, int qBin, float& sigmay, float& sigmax, float& sy1, float& sy2, float& sx1, float& sx2)
2305 
2306 {
2307  // Interpolate for a new set of track angles
2308 
2309  // Local variables
2310  int i;
2311  int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx, index;
2312  float yratio, yxratio, xxratio;
2313  float acotb, cotb;
2314  float yrms, xrms;
2315  bool flip_y;
2316 
2317 
2318  // Find the index corresponding to id
2319 
2320  index = -1;
2321  for(i=0; i<(int)thePixelTemp.size(); ++i) {
2322 
2323  if(id == thePixelTemp[i].head.ID) {
2324 
2325  index = i;
2326  break;
2327  }
2328  }
2329 
2330 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2331  if(index < 0 || index >= (int)thePixelTemp.size()) {
2332  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors can't find needed template ID = " << id << std::endl;
2333  }
2334 #else
2335  assert(index >= 0 && index < (int)thePixelTemp.size());
2336 #endif
2337 
2338 
2339 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2340  if(qBin < 0 || qBin > 5) {
2341  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors called with illegal qBin = " << qBin << std::endl;
2342  }
2343 #else
2344  assert(qBin >= 0 && qBin < 6);
2345 #endif
2346 
2347 // The error information for qBin > 3 is taken to be the same as qBin=3
2348 
2349  if(qBin > 3) {qBin = 3;}
2350  //
2351 
2352  // Interpolate the absolute value of cot(beta)
2353 
2354  acotb = fabs((double)cotbeta);
2355  cotb = cotbeta;
2356 
2357  // for some cosmics, the ususal gymnastics are incorrect
2358 
2359 // if(thePixelTemp[index].head.Dtype == 0) {
2360  cotb = acotb;
2361  flip_y = false;
2362  if(cotbeta < 0.) {flip_y = true;}
2363 // } else {
2364 // if(locBz < 0.) {
2365 // cotb = cotbeta;
2366 // flip_y = false;
2367 // } else {
2368 // cotb = -cotbeta;
2369 // flip_y = true;
2370 // }
2371 // }
2372 
2373  // Copy the charge scaling factor to the private variable
2374 
2375  Ny = thePixelTemp[index].head.NTy;
2376  Nyx = thePixelTemp[index].head.NTyx;
2377  Nxx = thePixelTemp[index].head.NTxx;
2378 
2379 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2380  if(Ny < 2 || Nyx < 1 || Nxx < 2) {
2381  throw cms::Exception("DataCorrupt") << "template ID = " << id_current << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
2382  }
2383 #else
2384  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
2385 #endif
2386  imaxx = Nyx - 1;
2387  imidy = Nxx/2;
2388 
2389  // next, loop over all y-angle entries
2390 
2391  ilow = 0;
2392  yratio = 0.;
2393 
2394  if(cotb >= thePixelTemp[index].enty[Ny-1].cotbeta) {
2395 
2396  ilow = Ny-2;
2397  yratio = 1.;
2398 
2399  } else {
2400 
2401  if(cotb >= thePixelTemp[index].enty[0].cotbeta) {
2402 
2403  for (i=0; i<Ny-1; ++i) {
2404 
2405  if( thePixelTemp[index].enty[i].cotbeta <= cotb && cotb < thePixelTemp[index].enty[i+1].cotbeta) {
2406 
2407  ilow = i;
2408  yratio = (cotb - thePixelTemp[index].enty[i].cotbeta)/(thePixelTemp[index].enty[i+1].cotbeta - thePixelTemp[index].enty[i].cotbeta);
2409  break;
2410  }
2411  }
2412  }
2413  }
2414 
2415  ihigh=ilow + 1;
2416 
2417  // Interpolate/store all y-related quantities (flip displacements when flip_y)
2418 
2419  sy1 = (1. - yratio)*thePixelTemp[index].enty[ilow].syone + yratio*thePixelTemp[index].enty[ihigh].syone;
2420  sy2 = (1. - yratio)*thePixelTemp[index].enty[ilow].sytwo + yratio*thePixelTemp[index].enty[ihigh].sytwo;
2421  yrms=(1. - yratio)*thePixelTemp[index].enty[ilow].yrms[qBin] + yratio*thePixelTemp[index].enty[ihigh].yrms[qBin];
2422 
2423 
2424  // next, loop over all x-angle entries, first, find relevant y-slices
2425 
2426  iylow = 0;
2427  yxratio = 0.;
2428 
2429  if(acotb >= thePixelTemp[index].entx[Nyx-1][0].cotbeta) {
2430 
2431  iylow = Nyx-2;
2432  yxratio = 1.;
2433 
2434  } else if(acotb >= thePixelTemp[index].entx[0][0].cotbeta) {
2435 
2436  for (i=0; i<Nyx-1; ++i) {
2437 
2438  if( thePixelTemp[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp[index].entx[i+1][0].cotbeta) {
2439 
2440  iylow = i;
2441  yxratio = (acotb - thePixelTemp[index].entx[i][0].cotbeta)/(thePixelTemp[index].entx[i+1][0].cotbeta - thePixelTemp[index].entx[i][0].cotbeta);
2442  break;
2443  }
2444  }
2445  }
2446 
2447  iyhigh=iylow + 1;
2448 
2449  ilow = 0;
2450  xxratio = 0.;
2451 
2452  if(cotalpha >= thePixelTemp[index].entx[0][Nxx-1].cotalpha) {
2453 
2454  ilow = Nxx-2;
2455  xxratio = 1.;
2456 
2457  } else {
2458 
2459  if(cotalpha >= thePixelTemp[index].entx[0][0].cotalpha) {
2460 
2461  for (i=0; i<Nxx-1; ++i) {
2462 
2463  if( thePixelTemp[index].entx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp[index].entx[0][i+1].cotalpha) {
2464 
2465  ilow = i;
2466  xxratio = (cotalpha - thePixelTemp[index].entx[0][i].cotalpha)/(thePixelTemp[index].entx[0][i+1].cotalpha - thePixelTemp[index].entx[0][i].cotalpha);
2467  break;
2468  }
2469  }
2470  }
2471  }
2472 
2473  ihigh=ilow + 1;
2474 
2475  sx1 = (1. - xxratio)*thePixelTemp[index].entx[0][ilow].sxone + xxratio*thePixelTemp[index].entx[0][ihigh].sxone;
2476  sx2 = (1. - xxratio)*thePixelTemp[index].entx[0][ilow].sxtwo + xxratio*thePixelTemp[index].entx[0][ihigh].sxtwo;
2477 
2478  xrms=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index].entx[iylow][ilow].xrms[qBin] + xxratio*thePixelTemp[index].entx[iylow][ihigh].xrms[qBin])
2479  +yxratio*((1. - xxratio)*thePixelTemp[index].entx[iyhigh][ilow].xrms[qBin] + xxratio*thePixelTemp[index].entx[iyhigh][ihigh].xrms[qBin]);
2480 
2481 
2482 
2483 
2484  // Take the errors and bias from the correct charge bin
2485 
2486  sigmay = yrms;
2487 
2488  sigmax = xrms;
2489 
2490  return;
2491 
2492 } // temperrors
2493 
2494 // ************************************************************************************************************
2504 // ************************************************************************************************************
2505 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)
2506 
2507 {
2508  // Interpolate for a new set of track angles
2509 
2510  // Local variables
2511  int i;
2512  int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx, index;
2513  float yratio, yxratio, xxratio;
2514  float acotb, cotb;
2515  float qfrac[4];
2516  bool flip_y;
2517 
2518  // Find the index corresponding to id
2519 
2520  index = -1;
2521  for(i=0; i<(int)thePixelTemp.size(); ++i) {
2522 
2523  if(id == thePixelTemp[i].head.ID) {
2524 
2525  index = i;
2526 // id_current = id;
2527  break;
2528  }
2529  }
2530 
2531 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2532  if(index < 0 || index >= (int)thePixelTemp.size()) {
2533  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::temperrors can't find needed template ID = " << id << std::endl;
2534  }
2535 #else
2536  assert(index >= 0 && index < (int)thePixelTemp.size());
2537 #endif
2538 
2539  //
2540 
2541  // Interpolate the absolute value of cot(beta)
2542 
2543  acotb = fabs((double)cotbeta);
2544  cotb = cotbeta;
2545 
2546 
2547  // for some cosmics, the ususal gymnastics are incorrect
2548 
2549 // if(thePixelTemp[index].head.Dtype == 0) {
2550  cotb = acotb;
2551  flip_y = false;
2552  if(cotbeta < 0.) {flip_y = true;}
2553 // } else {
2554 // if(locBz < 0.) {
2555 // cotb = cotbeta;
2556 // flip_y = false;
2557 // } else {
2558 // cotb = -cotbeta;
2559 // flip_y = true;
2560 // }
2561 // }
2562 
2563  // Copy the charge scaling factor to the private variable
2564 
2565  Ny = thePixelTemp[index].head.NTy;
2566  Nyx = thePixelTemp[index].head.NTyx;
2567  Nxx = thePixelTemp[index].head.NTxx;
2568 
2569 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2570  if(Ny < 2 || Nyx < 1 || Nxx < 2) {
2571  throw cms::Exception("DataCorrupt") << "template ID = " << id_current << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
2572  }
2573 #else
2574  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
2575 #endif
2576  imaxx = Nyx - 1;
2577  imidy = Nxx/2;
2578 
2579  // next, loop over all y-angle entries
2580 
2581  ilow = 0;
2582  yratio = 0.;
2583 
2584  if(cotb >= thePixelTemp[index].enty[Ny-1].cotbeta) {
2585 
2586  ilow = Ny-2;
2587  yratio = 1.;
2588 
2589  } else {
2590 
2591  if(cotb >= thePixelTemp[index].enty[0].cotbeta) {
2592 
2593  for (i=0; i<Ny-1; ++i) {
2594 
2595  if( thePixelTemp[index].enty[i].cotbeta <= cotb && cotb < thePixelTemp[index].enty[i+1].cotbeta) {
2596 
2597  ilow = i;
2598  yratio = (cotb - thePixelTemp[index].enty[i].cotbeta)/(thePixelTemp[index].enty[i+1].cotbeta - thePixelTemp[index].enty[i].cotbeta);
2599  break;
2600  }
2601  }
2602  }
2603  }
2604 
2605  ihigh=ilow + 1;
2606 
2607  // Interpolate/store all y-related quantities (flip displacements when flip_y)
2608  ny1_frac = (1. - yratio)*thePixelTemp[index].enty[ilow].fracyone + yratio*thePixelTemp[index].enty[ihigh].fracyone;
2609  ny2_frac = (1. - yratio)*thePixelTemp[index].enty[ilow].fracytwo + yratio*thePixelTemp[index].enty[ihigh].fracytwo;
2610 
2611  // next, loop over all x-angle entries, first, find relevant y-slices
2612 
2613  iylow = 0;
2614  yxratio = 0.;
2615 
2616  if(acotb >= thePixelTemp[index].entx[Nyx-1][0].cotbeta) {
2617 
2618  iylow = Nyx-2;
2619  yxratio = 1.;
2620 
2621  } else if(acotb >= thePixelTemp[index].entx[0][0].cotbeta) {
2622 
2623  for (i=0; i<Nyx-1; ++i) {
2624 
2625  if( thePixelTemp[index].entx[i][0].cotbeta <= acotb && acotb < thePixelTemp[index].entx[i+1][0].cotbeta) {
2626 
2627  iylow = i;
2628  yxratio = (acotb - thePixelTemp[index].entx[i][0].cotbeta)/(thePixelTemp[index].entx[i+1][0].cotbeta - thePixelTemp[index].entx[i][0].cotbeta);
2629  break;
2630  }
2631  }
2632  }
2633 
2634  iyhigh=iylow + 1;
2635 
2636  ilow = 0;
2637  xxratio = 0.;
2638 
2639  if(cotalpha >= thePixelTemp[index].entx[0][Nxx-1].cotalpha) {
2640 
2641  ilow = Nxx-2;
2642  xxratio = 1.;
2643 
2644  } else {
2645 
2646  if(cotalpha >= thePixelTemp[index].entx[0][0].cotalpha) {
2647 
2648  for (i=0; i<Nxx-1; ++i) {
2649 
2650  if( thePixelTemp[index].entx[0][i].cotalpha <= cotalpha && cotalpha < thePixelTemp[index].entx[0][i+1].cotalpha) {
2651 
2652  ilow = i;
2653  xxratio = (cotalpha - thePixelTemp[index].entx[0][i].cotalpha)/(thePixelTemp[index].entx[0][i+1].cotalpha - thePixelTemp[index].entx[0][i].cotalpha);
2654  break;
2655  }
2656  }
2657  }
2658  }
2659 
2660  ihigh=ilow + 1;
2661 
2662  for(i=0; i<3; ++i) {
2663  qfrac[i]=(1. - yxratio)*((1. - xxratio)*thePixelTemp[index].entx[iylow][ilow].qbfrac[i] + xxratio*thePixelTemp[index].entx[iylow][ihigh].qbfrac[i])
2664  +yxratio*((1. - xxratio)*thePixelTemp[index].entx[iyhigh][ilow].qbfrac[i] + xxratio*thePixelTemp[index].entx[iyhigh][ihigh].qbfrac[i]);
2665  }
2666  nx1_frac = (1. - yxratio)*((1. - xxratio)*thePixelTemp[index].entx[iylow][ilow].fracxone + xxratio*thePixelTemp[index].entx[iylow][ihigh].fracxone)
2667  +yxratio*((1. - xxratio)*thePixelTemp[index].entx[iyhigh][ilow].fracxone + xxratio*thePixelTemp[index].entx[iyhigh][ihigh].fracxone);
2668  nx2_frac = (1. - yxratio)*((1. - xxratio)*thePixelTemp[index].entx[iylow][ilow].fracxtwo + xxratio*thePixelTemp[index].entx[iylow][ihigh].fracxtwo)
2669  +yxratio*((1. - xxratio)*thePixelTemp[index].entx[iyhigh][ilow].fracxtwo + xxratio*thePixelTemp[index].entx[iyhigh][ihigh].fracxtwo);
2670 
2671 
2672 
2673  qbin_frac[0] = qfrac[0];
2674  qbin_frac[1] = qbin_frac[0] + qfrac[1];
2675  qbin_frac[2] = qbin_frac[1] + qfrac[2];
2676  qbin_frac[3] = 1.;
2677  return;
2678 
2679 } // qbin
2680 
2681 // *************************************************************************************************************************************
2683 
2689 // *************************************************************************************************************************************
2690 
2691 bool SiPixelTemplate::simpletemplate2D(float xhit, float yhit, std::vector<bool>& ydouble, std::vector<bool>& xdouble, float template2d[BXM2][BYM2])
2692 {
2693 
2694  // Local variables
2695 
2696  float x0, y0, xf, yf, xi, yi, sf, si, s0, qpix, slopey, slopex, ds;
2697  int i, j, jpix0, ipix0, jpixf, ipixf, jpix, ipix, nx, ny, anx, any, jmax, imax;
2698  float qtotal;
2699  // double path;
2700  std::list<SimplePixel> list;
2701  std::list<SimplePixel>::iterator listIter, listEnd;
2702 
2703  // Calculate the entry and exit points for the line charge from the track
2704 
2705  x0 = xhit - 0.5*pzsize*cota_current;
2706  y0 = yhit - 0.5*pzsize*cotb_current;
2707 
2708  jpix0 = floor(x0/pxsize)+1;
2709  ipix0 = floor(y0/pysize)+1;
2710 
2711  if(jpix0 < 0 || jpix0 > BXM3) {return false;}
2712  if(ipix0 < 0 || ipix0 > BYM3) {return false;}
2713 
2714  xf = xhit + 0.5*pzsize*cota_current + plorxwidth;
2715  yf = yhit + 0.5*pzsize*cotb_current + plorywidth;
2716 
2717  jpixf = floor(xf/pxsize)+1;
2718  ipixf = floor(yf/pysize)+1;
2719 
2720  if(jpixf < 0 || jpixf > BXM3) {return false;}
2721  if(ipixf < 0 || ipixf > BYM3) {return false;}
2722 
2723 // total charge length
2724 
2725  sf = sqrt((xf-x0)*(xf-x0) + (yf-y0)*(yf-y0));
2726  if((xf-x0) != 0.) {slopey = (yf-y0)/(xf-x0);} else { slopey = 1.e10;}
2727  if((yf-y0) != 0.) {slopex = (xf-x0)/(yf-y0);} else { slopex = 1.e10;}
2728 
2729 // use average charge in this direction
2730 
2731  qtotal = pqavg_avg;
2732 
2733  SimplePixel element;
2734  element.s = sf;
2735  element.x = xf;
2736  element.y = yf;
2737  element.i = ipixf;
2738  element.j = jpixf;
2739  element.btype = 0;
2740  list.push_back(element);
2741 
2742  // nx is the number of x interfaces crossed by the line charge
2743 
2744  nx = jpixf - jpix0;
2745  anx = abs(nx);
2746  if(anx > 0) {
2747  if(nx > 0) {
2748  for(j=jpix0; j<jpixf; ++j) {
2749  xi = pxsize*j;
2750  yi = slopey*(xi-x0) + y0;
2751  ipix = (int)(yi/pysize)+1;
2752  si = sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
2753  element.s = si;
2754  element.x = xi;
2755  element.y = yi;
2756  element.i = ipix;
2757  element.j = j;
2758  element.btype = 1;
2759  list.push_back(element);
2760  }
2761  } else {
2762  for(j=jpix0; j>jpixf; --j) {
2763  xi = pxsize*(j-1);
2764  yi = slopey*(xi-x0) + y0;
2765  ipix = (int)(yi/pysize)+1;
2766  si = sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
2767  element.s = si;
2768  element.x = xi;
2769  element.y = yi;
2770  element.i = ipix;
2771  element.j = j;
2772  element.btype = 1;
2773  list.push_back(element);
2774  }
2775  }
2776  }
2777 
2778  ny = ipixf - ipix0;
2779  any = abs(ny);
2780  if(any > 0) {
2781  if(ny > 0) {
2782  for(i=ipix0; i<ipixf; ++i) {
2783  yi = pysize*i;
2784  xi = slopex*(yi-y0) + x0;
2785  jpix = (int)(xi/pxsize)+1;
2786  si = sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
2787  element.s = si;
2788  element.x = xi;
2789  element.y = yi;
2790  element.i = i;
2791  element.j = jpix;
2792  element.btype = 2;
2793  list.push_back(element);
2794  }
2795  } else {
2796  for(i=ipix0; i>ipixf; --i) {
2797  yi = pysize*(i-1);
2798  xi = slopex*(yi-y0) + x0;
2799  jpix = (int)(xi/pxsize)+1;
2800  si = sqrt((xi-x0)*(xi-x0) + (yi-y0)*(yi-y0));
2801  element.s = si;
2802  element.x = xi;
2803  element.y = yi;
2804  element.i = i;
2805  element.j = jpix;
2806  element.btype = 2;
2807  list.push_back(element);
2808  }
2809  }
2810  }
2811 
2812  imax = std::max(ipix0, ipixf);
2813  jmax = std::max(jpix0, jpixf);
2814 
2815  // Sort the list according to the distance from the initial point
2816 
2817  list.sort();
2818 
2819  // Look for double pixels and adjust the list appropriately
2820 
2821  for(i=1; i<imax; ++i) {
2822  if(ydouble[i-1]) {
2823  listIter = list.begin();
2824  if(ny > 0) {
2825  while(listIter != list.end()) {
2826  if(listIter->i == i && listIter->btype == 2) {
2827  listIter = list.erase(listIter);
2828  continue;
2829  }
2830  if(listIter->i > i) {
2831  --(listIter->i);
2832  }
2833  ++listIter;
2834  }
2835  } else {
2836  while(listIter != list.end()) {
2837  if(listIter->i == i+1 && listIter->btype == 2) {
2838  listIter = list.erase(listIter);
2839  continue;
2840  }
2841  if(listIter->i > i+1) {
2842  --(listIter->i);
2843  }
2844  ++listIter;
2845  }
2846  }
2847  }
2848  }
2849 
2850  for(j=1; j<jmax; ++j) {
2851  if(xdouble[j-1]) {
2852  listIter = list.begin();
2853  if(nx > 0) {
2854  while(listIter != list.end()) {
2855  if(listIter->j == j && listIter->btype == 1) {
2856  listIter = list.erase(listIter);
2857  continue;
2858  }
2859  if(listIter->j > j) {
2860  --(listIter->j);
2861  }
2862  ++listIter;
2863  }
2864  } else {
2865  while(listIter != list.end()) {
2866  if(listIter->j == j+1 && listIter->btype == 1) {
2867  listIter = list.erase(listIter);
2868  continue;
2869  }
2870  if(listIter->j > j+1) {
2871  --(listIter->j);
2872  }
2873  ++listIter;
2874  }
2875  }
2876  }
2877  }
2878 
2879  // 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.
2880 
2881  s0 = 0.;
2882  listIter = list.begin();
2883  listEnd = list.end();
2884  for( ;listIter != listEnd; ++listIter) {
2885  si = listIter->s;
2886  ds = si - s0;
2887  s0 = si;
2888  j = listIter->j;
2889  i = listIter->i;
2890  if(sf > 0.) { qpix = qtotal*ds/sf;} else {qpix = qtotal;}
2891  template2d[j][i] += qpix;
2892  }
2893 
2894  return true;
2895 
2896 } // simpletemplate2D
2897 
2898 
2899 // ************************************************************************************************************
2904 // ************************************************************************************************************
2905 void SiPixelTemplate::vavilov_pars(double& mpv, double& sigma, double& kappa)
2906 
2907 {
2908  // Local variables
2909  int i;
2910  int ilow, ihigh, Ny;
2911  float yratio, cotb;
2912 
2913 // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta)
2914 
2915  cotb = sqrt(cotb_current*cotb_current + cota_current*cota_current);
2916 
2917 // Copy the charge scaling factor to the private variable
2918 
2919  Ny = thePixelTemp[index_id].head.NTy;
2920 
2921 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
2922  if(Ny < 2) {
2923  throw cms::Exception("DataCorrupt") << "template ID = " << id_current << "has too few entries: Ny = " << Ny << std::endl;
2924  }
2925 #else
2926  assert(Ny > 1);
2927 #endif
2928 
2929 // next, loop over all y-angle entries
2930 
2931  ilow = 0;
2932  yratio = 0.;
2933 
2934  if(cotb >= thePixelTemp[index_id].enty[Ny-1].cotbeta) {
2935 
2936  ilow = Ny-2;
2937  yratio = 1.;
2938 
2939  } else {
2940 
2941  if(cotb >= thePixelTemp[index_id].enty[0].cotbeta) {
2942 
2943  for (i=0; i<Ny-1; ++i) {
2944 
2945  if( thePixelTemp[index_id].enty[i].cotbeta <= cotb && cotb < thePixelTemp[index_id].enty[i+1].cotbeta) {
2946 
2947  ilow = i;
2948  yratio = (cotb - thePixelTemp[index_id].enty[i].cotbeta)/(thePixelTemp[index_id].enty[i+1].cotbeta - thePixelTemp[index_id].enty[i].cotbeta);
2949  break;
2950  }
2951  }
2952  }
2953  }
2954 
2955  ihigh=ilow + 1;
2956 
2957 // Interpolate Vavilov parameters
2958 
2959  pmpvvav = (1. - yratio)*thePixelTemp[index_id].enty[ilow].mpvvav + yratio*thePixelTemp[index_id].enty[ihigh].mpvvav;
2960  psigmavav = (1. - yratio)*thePixelTemp[index_id].enty[ilow].sigmavav + yratio*thePixelTemp[index_id].enty[ihigh].sigmavav;
2961  pkappavav = (1. - yratio)*thePixelTemp[index_id].enty[ilow].kappavav + yratio*thePixelTemp[index_id].enty[ihigh].kappavav;
2962 
2963 // Copy to parameter list
2964 
2965 
2966  mpv = (double)pmpvvav;
2967  sigma = (double)psigmavav;
2968  kappa = (double)pkappavav;
2969 
2970  return;
2971 
2972 } // vavilov_pars
2973 
SiPixelTemplateEntry::qavg_avg
float qavg_avg
average cluster charge of clusters that are less than qavg (normalize 2-D simple templates) ...
Definition: SiPixelTemplate.h:150
SiPixelTemplateEntry::runnum
int runnum
&lt; Basic template entry corresponding to a single set of track angles
Definition: SiPixelTemplate.h:82
SiPixelTemplateEntry::xflpar
float xflpar[4][6]
Aqfl-parameterized x-correction in 4 charge bins.
Definition: SiPixelTemplate.h:120
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:117
RecoTauCommonJetSelections_cfi.nbins
tuple nbins
Definition: RecoTauCommonJetSelections_cfi.py:118
BXSIZE
#define BXSIZE
Definition: SiPixelTemplateDefs.h:26
LOGERROR
#define LOGERROR(x)
Definition: SiPixelTemplate.cc:84
SiPixelTemplateEntry::clslenx
float clslenx
cluster x-length in pixels at signal height sxmax/2
Definition: SiPixelTemplate.h:103
SiPixelTemplate::qbin
int qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float &pixmx, float &sigmay, float &deltay, float &sigmax, float &deltax, float &sy1, float &dy1, float &sy2, float &dy2, float &sx1, float &dx1, float &sx2, float &dx2, float &lorywidth, float &lorxwidth)
Definition: SiPixelTemplate.cc:1984
SiPixelTemplateEntry::xtemp
float xtemp[9][TXSIZE]
templates for x-reconstruction (binned over 1 central pixel)
Definition: SiPixelTemplate.h:110
SiPixelTemplateEntry::qavg
float qavg
average cluster charge for this set of track angles (now includes threshold effects) ...
Definition: SiPixelTemplate.h:88
SiPixelTemplate::ytemp3d
void ytemp3d(int nypix, array_3d &ytemplate)
Definition: SiPixelTemplate.cc:1779
SiPixelTemplateEntry::syone
float syone
rms for one pixel y-clusters
Definition: SiPixelTemplate.h:92
SiPixelTemplateEntry::chi2yavgone
float chi2yavgone
average y chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:125
SiPixelTemplateEntry::dyone
float dyone
mean offset/correction for one pixel y-clusters
Definition: SiPixelTemplate.h:91
BYSIZE
#define BYSIZE
Definition: SiPixelTemplateDefs.h:18
SiPixelTemplateEntry::sigmavav
float sigmavav
&quot;sigma&quot; scale fctor for Vavilov distribution
Definition: SiPixelTemplate.h:105
SiPixelTemplateEntry::fracxtwo
float fracxtwo
fraction of double pixel sample with xsize = 1
Definition: SiPixelTemplate.h:149
SiPixelTemplateEntry::yavggen
float yavggen[4]
generic algorithm: average y-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:137
SiPixelTemplateEntry::sxmax
float sxmax
average pixel signal for x-projection of cluster
Definition: SiPixelTemplate.h:93
SimplePixel::btype
int btype
type of boundary (0=end, 1 = x-boundary, 2 = y-boundary)
Definition: SimplePixel.h:22
SiPixelTemplateEntry::xavgc2m
float xavgc2m[4]
1st pass chi2 min search: average x-bias of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:133
TXSIZE
#define TXSIZE
Definition: SiPixelTemplateDefs.h:23
SiPixelTemplateEntry::chi2xavgone
float chi2xavgone
average x chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:127
SiPixelTemplateEntry::fracytwo
float fracytwo
fraction of double pixel sample with ysize = 1
Definition: SiPixelTemplate.h:148
SiPixelTemplateEntry::yrms
float yrms[4]
average y-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:112
SiPixelTemplateHeader::NTxx
int NTxx
number of Template x-entries in each slice
Definition: SiPixelTemplate.h:165
BXM3
#define BXM3
Definition: SiPixelTemplateDefs.h:30
SiPixelTemplateEntry::ytemp
float ytemp[9][TYSIZE]
templates for y-reconstruction (binned over 1 central pixel)
Definition: SiPixelTemplate.h:108
BXM1
#define BXM1
Definition: SiPixelTemplateDefs.h:28
SiPixelTemplateEntry::ygx0c2m
float ygx0c2m[4]
1st pass chi2 min search: average y0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:131
SiPixelTemplateEntry::ygx0gen
float ygx0gen[4]
generic algorithm: average y0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:139
SiPixelTemplateEntry::xavg
float xavg[4]
average x-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:116
EcalCondDB.db
tuple db
Definition: EcalCondDB.py:151
abs
#define abs(x)
Definition: mlp_lapack.h:159
SimplePixel::y
float y
y coordinate of boundary intersection
Definition: SimplePixel.h:19
SiPixelTemplateHeader::fluence
float fluence
radiation fluence in n_eq/cm^2
Definition: SiPixelTemplate.h:169
SiPixelTemplateEntry::dytwo
float dytwo
mean offset/correction for one double-pixel y-clusters
Definition: SiPixelTemplate.h:96
SiPixelTemplateHeader::ysize
float ysize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate.h:175
SiPixelTemplateEntry::dxone
float dxone
mean offset/correction for one pixel x-clusters
Definition: SiPixelTemplate.h:94
SiPixelTemplateEntry::pixmax
float pixmax
maximum charge for individual pixels in cluster
Definition: SiPixelTemplate.h:89
SiPixelTemplateHeader::qscale
float qscale
Charge scaling to match cmssw and pixelav.
Definition: SiPixelTemplate.h:170
SiPixelTemplateEntry::qmin
float qmin
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
Definition: SiPixelTemplate.h:100
SiPixelTemplate::xflcorr
float xflcorr(int binq, float qflx)
Definition: SiPixelTemplate.cc:1573
SiPixelTemplateEntry::xgsiggen
float xgsiggen[4]
generic algorithm: average sigma_x from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:144
SiPixelTemplateHeader::ID
int ID
template ID number
Definition: SiPixelTemplate.h:160
SiPixelTemplateHeader::Bfield
float Bfield
Bfield in Tesla.
Definition: SiPixelTemplate.h:162
SiPixelTemplateEntry::cotalpha
float cotalpha
cot(alpha) is proportional to cluster length in x and is basis of interpolation
Definition: SiPixelTemplate.h:84
SiPixelTemplateHeader::temperature
float temperature
detector temperature in deg K
Definition: SiPixelTemplate.h:168
SiPixelTemplate::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBz)
Definition: SiPixelTemplate.cc:844
SiPixelTemplateEntry::yrmsc2m
float yrmsc2m[4]
1st pass chi2 min search: average y-rms of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:130
SiPixelTemplateHeader::lorywidth
float lorywidth
estimate of y-lorentz width from single pixel offset
Definition: SiPixelTemplate.h:172
BYM1
#define BYM1
Definition: SiPixelTemplateDefs.h:20
BXM2
#define BXM2
Definition: SiPixelTemplateDefs.h:29
SiPixelTemplateReco::array_2d
boost::multi_array< float, 2 > array_2d
Definition: SiPixelTemplateReco.h:69
SiPixelTemplateHeader::templ_version
int templ_version
Version number of the template to ensure code compatibility.
Definition: SiPixelTemplate.h:161
SiPixelTemplateEntry::xpar
float xpar[2][5]
projected x-pixel uncertainty parameterization
Definition: SiPixelTemplate.h:109
SiPixelTemplateEntry::ygx0
float ygx0[4]
average y0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:113
SiPixelTemplate::simpletemplate2D
bool simpletemplate2D(float xhitp, float yhitp, std::vector< bool > &ydouble, std::vector< bool > &xdouble, float template2d[BXM2][BYM2])
Make simple 2-D templates from track angles set in interpolate and hit position.
Definition: SiPixelTemplate.cc:2691
SiPixelTemplateEntry::xgx0gen
float xgx0gen[4]
generic algorithm: average x0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:143
SiPixelTemplate::xtemp3d
void xtemp3d(int nxpix, array_3d &xtemplate)
Definition: SiPixelTemplate.cc:1871
SiPixelTemplateEntry::symax
float symax
average pixel signal for y-projection of cluster
Definition: SiPixelTemplate.h:90
SiPixelTemplateEntry::xgsig
float xgsig[4]
average sigma_x from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:119
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
SiPixelTemplateEntry::cotbeta
float cotbeta
cot(beta) is proportional to cluster length in y and is basis of interpolation
Definition: SiPixelTemplate.h:86
SiPixelTemplateEntry::beta
float beta
beta track angle (defined in CMS CMS IN 2004/014)
Definition: SiPixelTemplate.h:85
SiPixelTemplate::vavilov_pars
void vavilov_pars(double &mpv, double &sigma, double &kappa)
Definition: SiPixelTemplate.cc:2905
SiPixelTemplate::xtemp
void xtemp(int fxbin, int lxbin, float xtemplate[41][BXSIZE])
Definition: SiPixelTemplate.cc:1703
ENDL
#define ENDL
Definition: SiPixelTemplate.cc:86
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
BHX
#define BHX
Definition: SiPixelTemplateDefs.h:27
SiPixelTemplate::yflcorr
float yflcorr(int binq, float qfly)
Definition: SiPixelTemplate.cc:1521
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:134
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:2304
SiPixelTemplateHeader::lorxwidth
float lorxwidth
estimate of x-lorentz width from single pixel offset
Definition: SiPixelTemplate.h:173
SiPixelTemplateEntry::xgx0
float xgx0[4]
average x0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:118
j
int j
Definition: DBlmapReader.cc:9
SiPixelTemplateEntry::ygsigc2m
float ygsigc2m[4]
1st pass chi2 min search: average sigma_y from Gaussian fit binned in 4 charge bins ...
Definition: SiPixelTemplate.h:132
SiPixelTemplateHeader::title
char title[80]
&lt; template header structure
Definition: SiPixelTemplate.h:159
BYM2
#define BYM2
Definition: SiPixelTemplateDefs.h:21
SimplePixel::j
int j
y index of traversed pixel
Definition: SimplePixel.h:21
f
double f[11][100]
Definition: MuScleFitUtils.cc:79
SiPixelTemplate::ytemp
void ytemp(int fybin, int lybin, float ytemplate[41][BYSIZE])
Definition: SiPixelTemplate.cc:1625
SiPixelTemplateEntry::ygsiggen
float ygsiggen[4]
generic algorithm: average sigma_y from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:140
SiPixelTemplateEntry::xgsigc2m
float xgsigc2m[4]
1st pass chi2 min search: average sigma_x from Gaussian fit binned in 4 charge bins ...
Definition: SiPixelTemplate.h:136
SiPixelTemplateEntry::fracxone
float fracxone
fraction of sample with xsize = 1
Definition: SiPixelTemplate.h:147
SiPixelTemplateEntry::mpvvav
float mpvvav
most probable charge in Vavilov distribution (not actually for larger kappa)
Definition: SiPixelTemplate.h:104
BYM3
#define BYM3
Definition: SiPixelTemplateDefs.h:22
SiPixelTemplateEntry::sxtwo
float sxtwo
rms for one double-pixel x-clusters
Definition: SiPixelTemplate.h:99
SiPixelTemplateHeader::Dtype
int Dtype
detector type (0=BPix, 1=FPix)
Definition: SiPixelTemplate.h:166
SiPixelTemplateDBObject::sVector
std::vector< float > sVector() const
Definition: SiPixelTemplateDBObject.h:65
SiPixelTemplateEntry::chi2yavg
float chi2yavg[4]
average y chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:121
SiPixelTemplateEntry::alpha
float alpha
alpha track angle (defined in CMS CMS IN 2004/014)
Definition: SiPixelTemplate.h:83
SiPixelTemplateStore::head
SiPixelTemplateHeader head
&lt; template storage structure
Definition: SiPixelTemplate.h:182
SiPixelTemplateEntry::chi2xminone
float chi2xminone
minimum of x chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:128
SiPixelTemplateEntry::ypar
float ypar[2][5]
projected y-pixel uncertainty parameterization
Definition: SiPixelTemplate.h:107
SiPixelTemplateEntry::qavg_spare
float qavg_spare
spare cluster charge
Definition: SiPixelTemplate.h:151
gen::k
int k[5][pyjets_maxn]
Definition: Pythia6Hadronizer.cc:86
SiPixelTemplate::pushfile
bool pushfile(int filenum)
Definition: SiPixelTemplate.cc:103
SiPixelTemplateHeader::s50
float s50
1/2 of the readout threshold in ADC units
Definition: SiPixelTemplate.h:171
TYSIZE
#define TYSIZE
Definition: SiPixelTemplateDefs.h:14
SiPixelTemplateHeader::zsize
float zsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate.h:176
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:2505
SiPixelTemplateHeader::Vbias
float Vbias
detector bias potential in Volts
Definition: SiPixelTemplate.h:167
SiPixelTemplateEntry::qbfrac
float qbfrac[3]
fraction of sample in qbin = 0-2 (&gt;=3 is the complement)
Definition: SiPixelTemplate.h:145
SiPixelTemplateHeader::xsize
float xsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate.h:174
SiPixelTemplateEntry::fracyone
float fracyone
fraction of sample with ysize = 1
Definition: SiPixelTemplate.h:146
SiPixelTemplateEntry::qmin2
float qmin2
tighter minimum cluster charge for valid hit (keeps 99.8% of simulated hits)
Definition: SiPixelTemplate.h:101
SiPixelTemplateEntry::xgx0c2m
float xgx0c2m[4]
1st pass chi2 min search: average x0 from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:135
SiPixelTemplateEntry::yavg
float yavg[4]
average y-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:111
BHY
#define BHY
Definition: SiPixelTemplateDefs.h:19
SiPixelTemplateEntry::chi2ymin
float chi2ymin[4]
minimum of y chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:122
SiPixelTemplateEntry::yavgc2m
float yavgc2m[4]
1st pass chi2 min search: average y-bias of reconstruction binned in 4 charge bins ...
Definition: SiPixelTemplate.h:129
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:163
SiPixelTemplate::xsigma2
void xsigma2(int fxpix, int lxpix, float sxthr, float xsum[BXSIZE], float xsig2[BXSIZE])
Definition: SiPixelTemplate.cc:1414
SiPixelTemplateEntry::dxtwo
float dxtwo
mean offset/correction for one double-pixel x-clusters
Definition: SiPixelTemplate.h:98
LOGINFO
#define LOGINFO(x)
Definition: SiPixelTemplate.cc:85
SimplePixel::i
int i
x index of traversed pixel
Definition: SimplePixel.h:20
SiPixelTemplate::ysigma2
void ysigma2(int fypix, int lypix, float sythr, float ysum[BYSIZE], float ysig2[BYSIZE])
SiPixelTemplateHeader::NTyx
int NTyx
number of Template y-slices of x entries
Definition: SiPixelTemplate.h:164
SiPixelTemplateEntry::yflpar
float yflpar[4][6]
Aqfl-parameterized y-correction in 4 charge bins.
Definition: SiPixelTemplate.h:115
SiPixelTemplateEntry::sxone
float sxone
rms for one pixel x-clusters
Definition: SiPixelTemplate.h:95
SiPixelTemplateEntry::chi2xmin
float chi2xmin[4]
minimum of x chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:124
SiPixelTemplateEntry::chi2yminone
float chi2yminone
minimum of y chi^2 for 1 pixel clusters
Definition: SiPixelTemplate.h:126
SiPixelTemplateEntry::costrk
float costrk[3]
direction cosines of tracks used to generate this entry
Definition: SiPixelTemplate.h:87
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:184
SiPixelTemplateEntry::ygsig
float ygsig[4]
average sigma_y from Gaussian fit binned in 4 charge bins
Definition: SiPixelTemplate.h:114
SiPixelTemplateStore::enty
SiPixelTemplateEntry enty[60]
60 Barrel y templates spanning cluster lengths from 0px to +18px [28 entries for fpix] ...
Definition: SiPixelTemplate.h:183
SiPixelTemplateEntry::sytwo
float sytwo
rms for one double-pixel y-clusters
Definition: SiPixelTemplate.h:97
SiPixelTemplateEntry::xrmsgen
float xrmsgen[4]
generic algorithm: average x-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:142
SiPixelTemplateEntry::chi2xavg
float chi2xavg[4]
average x chi^2 in 4 charge bins
Definition: SiPixelTemplate.h:123
array_3d
boost::multi_array< float, 3 > array_3d
Definition: SiPixelTemplate.h:79
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:106
SiPixelTemplateEntry::xavggen
float xavggen[4]
generic algorithm: average x-bias of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:141
SiPixelTemplateEntry::yrmsgen
float yrmsgen[4]
generic algorithm: average y-rms of reconstruction binned in 4 charge bins
Definition: SiPixelTemplate.h:138
SiPixelTemplateEntry::clsleny
float clsleny
cluster y-length in pixels at signal height symax/2
Definition: SiPixelTemplate.h:102
Generated for CMSSW Reference Manual by doxygen 1.8.5