CMS 3D CMS Logo
SiPixelTemplate2D.cc
Go to the documentation of this file.
1 //
2 // SiPixelTemplate2D.cc Version 2.65
3 //
4 // Full 2-D templates for cluster splitting, simulated cluster reweighting, and improved cluster probability
5 //
6 // Created by Morris Swartz on 12/01/09.
7 // 2009 __TheJohnsHopkinsUniversity__.
8 //
9 // V1.01 - fix qavg_ filling
10 // V1.02 - Add locBz to test if FPix use is out of range
11 // V1.03 - Fix edge checking on final template to increase template size and to properly truncate cluster
12 // v2.00 - Major changes to accommodate 2D reconstruction
13 // v2.10 - Change chi2 and error scaling information to work with partially reconstructed clusters
14 // v2.20 - Add cluster charge probability information, side loading for template generation
15 // v2.21 - Double derivative interval [improves fit convergence]
16 // v2.25 - Resize template store to accommodate FPix Templates
17 // v2.30 - Fix bug found by P. Shuetze that compromises sqlite file loading
18 // v2.35 - Add directory path selection to the ascii pushfile method
19 // v2.50 - Change template storage to dynamically allocated 2D arrays of SiPixelTemplateEntry2D structs
20 // v2.51 - Ensure that the derivative arrays are correctly zeroed between calls
21 // v2.52 - Improve cosmetics for increased style points from judges
22 // v2.60 - Fix FPix multiframe lookup problem [takes +-cotalpha and +-cotbeta]
23 // v2.61a - Code 2.60 fix correctly
24 // v2.65 - change double pixel flags to work with new shifted reco code definition
25 //
26 
27 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
28 #else
29 #include <math.h>
30 #endif
31 #include <algorithm>
32 #include <vector>
33 //#include "boost/multi_array.hpp"
34 #include <iostream>
35 #include <iomanip>
36 #include <sstream>
37 #include <fstream>
38 
39 
40 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
41 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplate2D.h"
42 #include "FWCore/ParameterSet/interface/FileInPath.h"
43 #include "FWCore/MessageLogger/interface/MessageLogger.h"
44 #include "FWCore/Utilities/interface/isFinite.h"
45 #define LOGERROR(x) LogError(x)
46 #define LOGINFO(x) LogInfo(x)
47 #define ENDL " "
48 #include "FWCore/Utilities/interface/Exception.h"
49 using namespace edm;
50 #else
51 #include "SiPixelTemplate2D.h"
52 #define LOGERROR(x) std::cout << x << ": "
53 #define LOGINFO(x) std::cout << x << ": "
54 #define ENDL std::endl
55 #endif
56 
57 //****************************************************************
62 //****************************************************************
63 bool SiPixelTemplate2D::pushfile(int filenum, std::vector< SiPixelTemplateStore2D > & pixelTemp, std::string dir)
64 {
65  // Add template stored in external file numbered filenum to theTemplateStore
66 
67  // Local variables
68  const int code_version={21};
69 
70  // Create a filename for this run
71  std::string tempfile = std::to_string(filenum);
72 
73 
74  // Create different path in CMSSW than standalone
75 
76 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
77  // If integer filenum has less than 4 digits, prepend 0's until we have four numerical characters, e.g. "0292"
78  int nzeros = 4-tempfile.length();
79  if (nzeros > 0)
80  tempfile = std::string(nzeros, '0') + tempfile;
82 
83  tempfile = dir + "template_summary2D_zp" + tempfile + ".out";
84  edm::FileInPath file( tempfile ); // Find the file in CMSSW
85  tempfile = file.fullPath(); // Put it back with the whole path.
86 
87 #else
88  // This is the same as above, but more elegant.
89  std::ostringstream tout;
90  tout << "template_summary2D_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
91  tempfile = tout.str();
92 
93 #endif
94 
95  // Open the template file
96  //
97  std::ifstream in_file(tempfile);
98  if(in_file.is_open() && in_file.good()) {
99 
100  // Create a local template storage entry
101  SiPixelTemplateStore2D theCurrentTemp;
102 
103  // Read-in a header string first and print it
104  char c;
105  for (int i=0; (c=in_file.get()) != '\n'; ++i) {
106  if(i < 79) {theCurrentTemp.head.title[i] = c;}
107  else {theCurrentTemp.head.title[79] ='\0';}
108  }
109  LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
110 
111  // next, the header information
112  in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
113  >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
114  >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
115 
116  if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 0A, no template load" << ENDL; return false;}
117 
118  if(theCurrentTemp.head.templ_version > 17) {
119  in_file >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias
120  >> theCurrentTemp.head.lorxbias >> theCurrentTemp.head.fbin[0]
121  >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
122 
123  if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 0B, no template load"
124  << ENDL; return false;}
125  } else {
126 // This is for older [legacy] payloads
127  theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
128  theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth/2.f;
129  theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth/2.f;
130  theCurrentTemp.head.fbin[0] = 1.5f;
131  theCurrentTemp.head.fbin[1] = 1.00f;
132  theCurrentTemp.head.fbin[2] = 0.85f;
133  }
134 
135  LOGINFO("SiPixelTemplate2D") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
136  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
137  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
138  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
139  << ", 1/2 multi dcol threshold " << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold " << theCurrentTemp.head.ss50
140  << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", y Lorentz Bias " << theCurrentTemp.head.lorybias
141  << ", x Lorentz width " << theCurrentTemp.head.lorxwidth << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
142  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", " << theCurrentTemp.head.fbin[1]
143  << ", " << theCurrentTemp.head.fbin[2]
144  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
145 
146  if(theCurrentTemp.head.templ_version < code_version) {LOGERROR("SiPixelTemplate2D") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
147 
148  if(theCurrentTemp.head.NTy != 0) {LOGERROR("SiPixelTemplate2D") << "Trying to load 1-d template info into the 2-d template object, check your DB/global tag!" << ENDL; return false;}
149 
150  // next, layout the 2-d structure needed to store template
151 
152  theCurrentTemp.entry = new SiPixelTemplateEntry2D*[theCurrentTemp.head.NTyx];
153  theCurrentTemp.entry[0] = new SiPixelTemplateEntry2D[theCurrentTemp.head.NTyx*theCurrentTemp.head.NTxx];
154  for(int i = 1; i < theCurrentTemp.head.NTyx; ++i) theCurrentTemp.entry[i] = theCurrentTemp.entry[i-1] + theCurrentTemp.head.NTxx;
155 
156  // Read in the file info
157 
158  for (int iy=0; iy < theCurrentTemp.head.NTyx; ++iy) {
159  for(int jx=0; jx < theCurrentTemp.head.NTxx; ++jx) {
160 
161  in_file >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0]
162  >> theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2];
163 
164  if(in_file.fail()) {
165  LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
166  delete [] theCurrentTemp.entry[0];
167  delete [] theCurrentTemp.entry;
168  return false;
169  }
170 
171  // Calculate cot(alpha) and cot(beta) for this entry
172 
173  theCurrentTemp.entry[iy][jx].cotalpha = theCurrentTemp.entry[iy][jx].costrk[0]/theCurrentTemp.entry[iy][jx].costrk[2];
174 
175  theCurrentTemp.entry[iy][jx].cotbeta = theCurrentTemp.entry[iy][jx].costrk[1]/theCurrentTemp.entry[iy][jx].costrk[2];
176 
177  in_file >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >> theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin
178  >> theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >> theCurrentTemp.entry[iy][jx].jxmax;
179 
180  if(in_file.fail()) {
181  LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
182  delete [] theCurrentTemp.entry[0];
183  delete [] theCurrentTemp.entry;
184  return false;
185  }
186 
187  for (int k=0; k<2; ++k) {
188 
189  in_file >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1]
190  >> theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >> theCurrentTemp.entry[iy][jx].xypar[k][4];
191 
192  if(in_file.fail()) {
193  LOGERROR("SiPixelTemplate2D") << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
194  delete [] theCurrentTemp.entry[0];
195  delete [] theCurrentTemp.entry;
196  return false;
197  }
198  }
199 
200  for (int k=0; k<2; ++k) {
201 
202  in_file >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1]
203  >> theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >> theCurrentTemp.entry[iy][jx].lanpar[k][4];
204 
205  if(in_file.fail()) {
206  LOGERROR("SiPixelTemplate2D") << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
207  delete [] theCurrentTemp.entry[0];
208  delete [] theCurrentTemp.entry;
209  return false;
210  }
211 
212  }
213 
214 
215 // Read the 2D template entries as floats [they are formatted that way] and cast to short ints
216 
217  float dummy[T2YSIZE];
218  for (int l=0; l<7; ++l) {
219  for (int k=0; k<7; ++k) {
220  for (int j=0; j<T2XSIZE; ++j) {
221  for (int i=0; i<T2YSIZE; ++i) {
222  in_file >> dummy[i];
223  }
224  if(in_file.fail()) {
225  LOGERROR("SiPixelTemplate2D") << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
226  delete [] theCurrentTemp.entry[0];
227  delete [] theCurrentTemp.entry;
228  return false;
229  }
230  for (int i=0; i<T2YSIZE; ++i) {
231  theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j] = (short int)dummy[i];
232  }
233  }
234  }
235  }
236 
237 
238  in_file >> theCurrentTemp.entry[iy][jx].chi2ppix >> theCurrentTemp.entry[iy][jx].chi2scale >> theCurrentTemp.entry[iy][jx].offsetx[0] >> theCurrentTemp.entry[iy][jx].offsetx[1]
239  >> theCurrentTemp.entry[iy][jx].offsetx[2] >> theCurrentTemp.entry[iy][jx].offsetx[3]>> theCurrentTemp.entry[iy][jx].offsety[0] >> theCurrentTemp.entry[iy][jx].offsety[1]
240  >> theCurrentTemp.entry[iy][jx].offsety[2] >> theCurrentTemp.entry[iy][jx].offsety[3];
241 
242  if(in_file.fail()) {
243  LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
244  delete [] theCurrentTemp.entry[0];
245  delete [] theCurrentTemp.entry;
246  return false;
247  }
248 
249  in_file >> theCurrentTemp.entry[iy][jx].clsleny >> theCurrentTemp.entry[iy][jx].clslenx >> theCurrentTemp.entry[iy][jx].mpvvav >> theCurrentTemp.entry[iy][jx].sigmavav
250  >> theCurrentTemp.entry[iy][jx].kappavav >> theCurrentTemp.entry[iy][jx].scalexavg >> theCurrentTemp.entry[iy][jx].scaleyavg >> theCurrentTemp.entry[iy][jx].delyavg >> theCurrentTemp.entry[iy][jx].delysig >> theCurrentTemp.entry[iy][jx].spare[0];
251 
252  if(in_file.fail()) {
253  LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
254  delete [] theCurrentTemp.entry[0];
255  delete [] theCurrentTemp.entry;
256  return false;
257  }
258 
259  in_file >> theCurrentTemp.entry[iy][jx].scalex[0] >> theCurrentTemp.entry[iy][jx].scalex[1] >> theCurrentTemp.entry[iy][jx].scalex[2] >> theCurrentTemp.entry[iy][jx].scalex[3]
260  >> theCurrentTemp.entry[iy][jx].scaley[0] >> theCurrentTemp.entry[iy][jx].scaley[1] >> theCurrentTemp.entry[iy][jx].scaley[2] >> theCurrentTemp.entry[iy][jx].scaley[3]
261  >> theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2];
262 
263  if(in_file.fail()) {
264  LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
265  delete [] theCurrentTemp.entry[0];
266  delete [] theCurrentTemp.entry;
267  return false;
268  }
269 
270  }
271 
272  }
273 
274 
275  in_file.close();
276 
277  // Add this template to the store
278 
279  pixelTemp.push_back(theCurrentTemp);
280 
281  return true;
282 
283  } else {
284 
285  // If file didn't open, report this
286 
287  LOGERROR("SiPixelTemplate2D") << "Error opening File" << tempfile << ENDL;
288  return false;
289 
290  }
291 
292 } // TempInit
293 
294 
295 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
296 
297 //****************************************************************
301 //****************************************************************
302 bool SiPixelTemplate2D::pushfile(const SiPixel2DTemplateDBObject& dbobject, std::vector< SiPixelTemplateStore2D > & pixelTemp)
303 {
304  // Add template stored in external dbobject to theTemplateStore
305 
306  const int code_version={21};
307 
308  // We must create a new object because dbobject must be a const and our stream must not be
309  SiPixel2DTemplateDBObject db = dbobject;
310 
311  // Create a local template storage entry
312  SiPixelTemplateStore2D theCurrentTemp;
313 
314  // Fill the template storage for each template calibration stored in the db
315  for(int m=0; m<db.numOfTempl(); ++m)
316  {
317 
318  // Read-in a header string first and print it
319 
320  SiPixel2DTemplateDBObject::char2float temp;
321  for (int i=0; i<20; ++i) {
322  temp.f = db.sVector()[db.index()];
323  theCurrentTemp.head.title[4*i] = temp.c[0];
324  theCurrentTemp.head.title[4*i+1] = temp.c[1];
325  theCurrentTemp.head.title[4*i+2] = temp.c[2];
326  theCurrentTemp.head.title[4*i+3] = temp.c[3];
327  db.incrementIndex(1);
328  }
329  theCurrentTemp.head.title[79] = '\0';
330  LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
331 
332  // next, the header information
333 
334  db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
335  >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
336  >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
337 
338  if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 0A, no template load" << ENDL; return false;}
339 
340  LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title
341  <<" code version = "<<code_version
342  <<" object version "<<theCurrentTemp.head.templ_version
343  << ENDL;
344 
345  if(theCurrentTemp.head.templ_version > 17) {
346  db >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >> theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
347 
348  if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 0B, no template load"
349  << ENDL; return false;}
350  } else {
351 // This is for older [legacy] payloads and the numbers are indeed magic [they are part of the payload for v>17]
352  theCurrentTemp.head.ss50 = theCurrentTemp.head.s50;
353  theCurrentTemp.head.lorybias = theCurrentTemp.head.lorywidth/2.f;
354  theCurrentTemp.head.lorxbias = theCurrentTemp.head.lorxwidth/2.f;
355  theCurrentTemp.head.fbin[0] = 1.50f;
356  theCurrentTemp.head.fbin[1] = 1.00f;
357  theCurrentTemp.head.fbin[2] = 0.85f;
358  }
359 
360  LOGINFO("SiPixelTemplate2D")
361  << "Template ID = " << theCurrentTemp.head.ID << ", Template Version "
362  << theCurrentTemp.head.templ_version << ", Bfield = "
363  << theCurrentTemp.head.Bfield<< ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = "
364  << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = "
365  << theCurrentTemp.head.Dtype<< ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
366  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence
367  << ", Q-scaling factor " << theCurrentTemp.head.qscale
368  << ", 1/2 multi dcol threshold " << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold "
369  << theCurrentTemp.head.ss50<< ", y Lorentz Width " << theCurrentTemp.head.lorywidth
370  << ", y Lorentz Bias " << theCurrentTemp.head.lorybias
371  << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
372  << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
373  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0]
374  << ", " << theCurrentTemp.head.fbin[1]
375  << ", " << theCurrentTemp.head.fbin[2]
376  << ", pixel x-size " << theCurrentTemp.head.xsize
377  << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
378 
379  if(theCurrentTemp.head.templ_version < code_version) {
380  LOGINFO("SiPixelTemplate2D") << "code expects version "<< code_version << " finds "
381  << theCurrentTemp.head.templ_version <<", load anyway " << ENDL;
382  }
383 
384  if(theCurrentTemp.head.NTy != 0) {LOGERROR("SiPixelTemplate2D") << "Trying to load 1-d template info into the 2-d template object, check your DB/global tag!" << ENDL; return false;}
385 
386 
387  // next, layout the 2-d structure needed to store template
388 
389  theCurrentTemp.entry = new SiPixelTemplateEntry2D*[theCurrentTemp.head.NTyx];
390  theCurrentTemp.entry[0] = new SiPixelTemplateEntry2D[theCurrentTemp.head.NTyx*theCurrentTemp.head.NTxx];
391  for(int i = 1; i < theCurrentTemp.head.NTyx; ++i) theCurrentTemp.entry[i] = theCurrentTemp.entry[i-1] + theCurrentTemp.head.NTxx;
392 
393  // Read in the file info
394 
395  for (int iy=0; iy < theCurrentTemp.head.NTyx; ++iy) {
396  for(int jx=0; jx < theCurrentTemp.head.NTxx; ++jx) {
397 
398  db >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0]
399  >> theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2];
400 
401  if(db.fail()) {
402  LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
403  delete [] theCurrentTemp.entry[0];
404  delete [] theCurrentTemp.entry;
405  return false;
406  }
407 
408  // Calculate cot(alpha) and cot(beta) for this entry
409 
410  theCurrentTemp.entry[iy][jx].cotalpha = theCurrentTemp.entry[iy][jx].costrk[0]/theCurrentTemp.entry[iy][jx].costrk[2];
411 
412  theCurrentTemp.entry[iy][jx].cotbeta = theCurrentTemp.entry[iy][jx].costrk[1]/theCurrentTemp.entry[iy][jx].costrk[2];
413 
414  db >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >> theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin
415  >> theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >> theCurrentTemp.entry[iy][jx].jxmax;
416 
417  if(db.fail()) {
418  LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
419  delete [] theCurrentTemp.entry[0];
420  delete [] theCurrentTemp.entry;
421  return false;
422  }
423 
424  for (int k=0; k<2; ++k) {
425 
426  db >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1]
427  >> theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >> theCurrentTemp.entry[iy][jx].xypar[k][4];
428 
429  if(db.fail()) {
430  LOGERROR("SiPixelTemplate2D") << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
431  delete [] theCurrentTemp.entry[0];
432  delete [] theCurrentTemp.entry;
433  return false;
434  }
435 
436  }
437 
438  for (int k=0; k<2; ++k) {
439 
440  db >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1]
441  >> theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >> theCurrentTemp.entry[iy][jx].lanpar[k][4];
442 
443  if(db.fail()) {
444  LOGERROR("SiPixelTemplate2D") << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
445  delete [] theCurrentTemp.entry[0];
446  delete [] theCurrentTemp.entry;
447  return false;
448  }
449 
450  }
451 
452 // Read the 2D template entries as floats [they are formatted that way] and cast to short ints
453 
454  float dummy[T2YSIZE];
455  for (int l=0; l<7; ++l) {
456  for (int k=0; k<7; ++k) {
457  for (int j=0; j<T2XSIZE; ++j) {
458  for (int i=0; i<T2YSIZE; ++i) {
459  db >> dummy[i];
460  }
461  if(db.fail()) {
462  LOGERROR("SiPixelTemplate2D") << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
463  delete [] theCurrentTemp.entry[0];
464  delete [] theCurrentTemp.entry;
465  return false;
466  }
467  for (int i=0; i<T2YSIZE; ++i) {
468  theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j] = (short int)dummy[i];
469  }
470  }
471  }
472  }
473 
474  db >> theCurrentTemp.entry[iy][jx].chi2ppix >> theCurrentTemp.entry[iy][jx].chi2scale >> theCurrentTemp.entry[iy][jx].offsetx[0] >> theCurrentTemp.entry[iy][jx].offsetx[1]
475  >> theCurrentTemp.entry[iy][jx].offsetx[2] >> theCurrentTemp.entry[iy][jx].offsetx[3]>> theCurrentTemp.entry[iy][jx].offsety[0] >> theCurrentTemp.entry[iy][jx].offsety[1]
476  >> theCurrentTemp.entry[iy][jx].offsety[2] >> theCurrentTemp.entry[iy][jx].offsety[3];
477 
478  if(db.fail()) {
479  LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
480  delete [] theCurrentTemp.entry[0];
481  delete [] theCurrentTemp.entry;
482  return false;
483  }
484 
485  db >> theCurrentTemp.entry[iy][jx].clsleny >> theCurrentTemp.entry[iy][jx].clslenx >> theCurrentTemp.entry[iy][jx].mpvvav >> theCurrentTemp.entry[iy][jx].sigmavav
486  >> theCurrentTemp.entry[iy][jx].kappavav >> theCurrentTemp.entry[iy][jx].scalexavg >> theCurrentTemp.entry[iy][jx].scaleyavg >> theCurrentTemp.entry[iy][jx].delyavg >> theCurrentTemp.entry[iy][jx].delysig >> theCurrentTemp.entry[iy][jx].spare[0];
487 
488  if(db.fail()) {
489  LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
490  delete [] theCurrentTemp.entry[0];
491  delete [] theCurrentTemp.entry;
492  return false;
493  }
494 
495  db >> theCurrentTemp.entry[iy][jx].scalex[0] >> theCurrentTemp.entry[iy][jx].scalex[1] >> theCurrentTemp.entry[iy][jx].scalex[2] >> theCurrentTemp.entry[iy][jx].scalex[3]
496  >> theCurrentTemp.entry[iy][jx].scaley[0] >> theCurrentTemp.entry[iy][jx].scaley[1] >> theCurrentTemp.entry[iy][jx].scaley[2] >> theCurrentTemp.entry[iy][jx].scaley[3]
497  >> theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2];
498 
499  if(db.fail()) {
500  LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL;
501  delete [] theCurrentTemp.entry[0];
502  delete [] theCurrentTemp.entry;
503  return false;
504  }
505 
506  }
507  }
508 
509 // Add this template to the store
510 
511  pixelTemp.push_back(theCurrentTemp);
512 
513  }
514 
515 
516  return true;
517 
518 } // TempInit
519 
520 #endif
521 
522 
523 bool SiPixelTemplate2D::getid(int id)
524 {
525  if(id != id_current_) {
526 
527  // Find the index corresponding to id
528 
529  index_id_ = -1;
530  for(int i=0; i<(int)thePixelTemp_.size(); ++i) {
531 
532  if(id == thePixelTemp_[i].head.ID) {
533 
534  index_id_ = i;
535  id_current_ = id;
536 
537  // Copy the charge scaling factor to the private variable
538 
539  Dtype_ = thePixelTemp_[index_id_].head.Dtype;
540 
541  // Copy the charge scaling factor to the private variable
542 
543  qscale_ = thePixelTemp_[index_id_].head.qscale;
544 
545  // Copy the pseudopixel signal size to the private variable
546 
547  s50_ = thePixelTemp_[index_id_].head.s50;
548 
549  // Copy Qbinning info to private variables
550 
551  for(int j=0; j<3; ++j) {fbin_[j] = thePixelTemp_[index_id_].head.fbin[j];}
552 
553  // Copy the Lorentz widths to private variables
554 
555  lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
556  lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
557 
558  // Copy the pixel sizes private variables
559 
560  xsize_ = thePixelTemp_[index_id_].head.xsize;
561  ysize_ = thePixelTemp_[index_id_].head.ysize;
562  zsize_ = thePixelTemp_[index_id_].head.zsize;
563 
564  // Determine the size of this template
565 
566  Nyx_ = thePixelTemp_[index_id_].head.NTyx;
567  Nxx_ = thePixelTemp_[index_id_].head.NTxx;
568 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
569  if(Nyx_ < 2 || Nxx_ < 2) {
570  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Nyx/Nxx = " << Nyx_ << "/" << Nxx_ << std::endl;
571  }
572 #else
573  assert(Nyx_ > 1 && Nxx_ > 1);
574 #endif
575  int imidx = Nxx_/2;
576 
577  cotalpha0_ = thePixelTemp_[index_id_].entry[0][0].cotalpha;
578  cotalpha1_ = thePixelTemp_[index_id_].entry[0][Nxx_-1].cotalpha;
579  deltacota_ = (cotalpha1_-cotalpha0_)/(float)(Nxx_-1);
580 
581  cotbeta0_ = thePixelTemp_[index_id_].entry[0][imidx].cotbeta;
582  cotbeta1_ = thePixelTemp_[index_id_].entry[Nyx_-1][imidx].cotbeta;
583  deltacotb_ = (cotbeta1_-cotbeta0_)/(float)(Nyx_-1);
584 
585  break;
586  }
587  }
588  }
589 
590 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
591  if(index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
592  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
593  << ", Are you using the correct global tag?" << std::endl;
594  }
595 #else
596  assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
597 #endif
598  return true;
599 }
600 
601 // *************************************************************************************************************************************
612 // *************************************************************************************************************************************
613 
614 bool SiPixelTemplate2D::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
615 {
616 
617 
618  // Interpolate for a new set of track angles
619 
620  // Local variables
621 
622  float acotb, dcota, dcotb;
623 
624  // Check to see if interpolation is valid
625 
626  if(id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
627  cota_current_ = cotalpha; cotb_current_ = cotbeta;
628  success_ = getid(id);
629  }
630 
631 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
632  if(index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
633  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
634  << ", Are you using the correct global tag?" << std::endl;
635  }
636 #else
637  assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
638 #endif
639 
640  // Check angle limits and et up interpolation parameters
641 
642  float cota = cotalpha;
643  flip_x_ = false;
644  flip_y_ = false;
645  switch(Dtype_) {
646  case 0:
647  if(cotbeta < 0.f) {flip_y_ = true;}
648  break;
649  case 1:
650  if(locBz > 0.f) {flip_y_ = true;}
651  break;
652  case 2:
653  case 3:
654  case 4:
655  case 5:
656  if(locBx*locBz < 0.f) {
657  cota = fabs(cotalpha);
658  flip_x_ = true;
659  }
660  if(locBx < 0.f) {flip_y_ = true;}
661  break;
662  default:
663 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
664  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
665 
666  //check for nan's
667  if (!edm::isFinite(cotalpha) || !edm::isFinite(cotbeta)) {
668  success_ = false;
669  return success_;
670  }
671 #else
672  std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
673 #endif
674  }
675 
676  if(cota < cotalpha0_) {
677  success_ = false;
678  jx0_ = 0;
679  jx1_ = 1;
680  adcota_ = 0.f;
681  } else if(cota > cotalpha1_) {
682  success_ = false;
683  jx0_ = Nxx_ - 1;
684  jx1_ = jx0_ - 1;
685  adcota_ = 0.f;
686  } else {
687  jx0_ = (int)((cota-cotalpha0_)/deltacota_+0.5f);
688  dcota = (cota - (cotalpha0_ + jx0_*deltacota_))/deltacota_;
689  adcota_ = fabs(dcota);
690  if(dcota > 0.f) {jx1_ = jx0_ + 1;if(jx1_ > Nxx_-1) jx1_ = jx0_-1;} else {jx1_ = jx0_ - 1; if(jx1_ < 0) jx1_ = jx0_+1;}
691  }
692 
693  // Interpolate the absolute value of cot(beta)
694 
695  acotb = fabs(cotbeta);
696 
697  if(acotb < cotbeta0_) {
698  success_ = false;
699  iy0_ = 0;
700  iy1_ = 1;
701  adcotb_ = 0.f;
702  } else if(acotb > cotbeta1_) {
703  success_ = false;
704  iy0_ = Nyx_ - 1;
705  iy1_ = iy0_ - 1;
706  adcotb_ = 0.f;
707  } else {
708  iy0_ = (int)((acotb-cotbeta0_)/deltacotb_+0.5f);
709  dcotb = (acotb - (cotbeta0_ + iy0_*deltacotb_))/deltacotb_;
710  adcotb_ = fabs(dcotb);
711  if(dcotb > 0.f) {iy1_ = iy0_ + 1; if(iy1_ > Nyx_-1) iy1_ = iy0_-1;} else {iy1_ = iy0_ - 1; if(iy1_ < 0) iy1_ = iy0_+1;}
712  }
713 
714  // Calculate signed quantities
715 
716  lorydrift_ = lorywidth_/2.;
717  if(flip_y_) lorydrift_ = -lorydrift_;
718  lorxdrift_ = lorxwidth_/2.;
719  if(flip_x_) lorxdrift_ = -lorxdrift_;
720 
721  // Use pointers to the three angle pairs used in the interpolation
722 
723 
724  entry00_ = &thePixelTemp_[index_id_].entry[iy0_][jx0_];
725  entry10_ = &thePixelTemp_[index_id_].entry[iy1_][jx0_];
726  entry01_ = &thePixelTemp_[index_id_].entry[iy0_][jx1_];
727 
728  // Interpolate things in cot(alpha)-cot(beta)
729 
730  qavg_ = entry00_->qavg
731  +adcota_*(entry01_->qavg - entry00_->qavg)
732  +adcotb_*(entry10_->qavg - entry00_->qavg);
733 
734  pixmax_ = entry00_->pixmax
735  +adcota_*(entry01_->pixmax - entry00_->pixmax)
736  +adcotb_*(entry10_->pixmax - entry00_->pixmax);
737 
738  sxymax_ = entry00_->sxymax
739  +adcota_*(entry01_->sxymax - entry00_->sxymax)
740  +adcotb_*(entry10_->sxymax - entry00_->sxymax);
741 
742  chi2avgone_ = entry00_->chi2avgone
743  +adcota_*(entry01_->chi2avgone - entry00_->chi2avgone)
744  +adcotb_*(entry10_->chi2avgone - entry00_->chi2avgone);
745 
746  chi2minone_ = entry00_->chi2minone
747  +adcota_*(entry01_->chi2minone - entry00_->chi2minone)
748  +adcotb_*(entry10_->chi2minone - entry00_->chi2minone);
749 
750  clsleny_ = entry00_->clsleny
751  +adcota_*(entry01_->clsleny - entry00_->clsleny)
752  +adcotb_*(entry10_->clsleny - entry00_->clsleny);
753 
754  clslenx_ = entry00_->clslenx
755  +adcota_*(entry01_->clslenx - entry00_->clslenx)
756  +adcotb_*(entry10_->clslenx - entry00_->clslenx);
757 
758 
759  chi2ppix_ = entry00_->chi2ppix
760  +adcota_*(entry01_->chi2ppix - entry00_->chi2ppix)
761  +adcotb_*(entry10_->chi2ppix - entry00_->chi2ppix);
762 
763  chi2scale_ = entry00_->chi2scale
764  +adcota_*(entry01_->chi2scale - entry00_->chi2scale)
765  +adcotb_*(entry10_->chi2scale - entry00_->chi2scale);
766 
767  scaleyavg_ = entry00_->scaleyavg
768  +adcota_*(entry01_->scaleyavg - entry00_->scaleyavg)
769  +adcotb_*(entry10_->scaleyavg - entry00_->scaleyavg);
770 
771  scalexavg_ = entry00_->scalexavg
772  +adcota_*(entry01_->scalexavg - entry00_->scalexavg)
773  +adcotb_*(entry10_->scalexavg - entry00_->scalexavg);
774 
775  delyavg_ = entry00_->delyavg
776  +adcota_*(entry01_->delyavg - entry00_->delyavg)
777  +adcotb_*(entry10_->delyavg - entry00_->delyavg);
778 
779  delysig_ = entry00_->delysig
780  +adcota_*(entry01_->delysig - entry00_->delysig)
781  +adcotb_*(entry10_->delysig - entry00_->delysig);
782 
783  mpvvav_ = entry00_->mpvvav
784  +adcota_*(entry01_->mpvvav - entry00_->mpvvav)
785  +adcotb_*(entry10_->mpvvav - entry00_->mpvvav);
786 
787  sigmavav_ = entry00_->sigmavav
788  +adcota_*(entry01_->sigmavav - entry00_->sigmavav)
789  +adcotb_*(entry10_->sigmavav - entry00_->sigmavav);
790 
791  kappavav_ = entry00_->kappavav
792  +adcota_*(entry01_->kappavav - entry00_->kappavav)
793  +adcotb_*(entry10_->kappavav - entry00_->kappavav);
794 
795  for(int i=0; i<4 ; ++i) {
796  scalex_[i] = entry00_->scalex[i]
797  +adcota_*(entry01_->scalex[i] - entry00_->scalex[i])
798  +adcotb_*(entry10_->scalex[i] - entry00_->scalex[i]);
799 
800  scaley_[i] = entry00_->scaley[i]
801  +adcota_*(entry01_->scaley[i] - entry00_->scaley[i])
802  +adcotb_*(entry10_->scaley[i] - entry00_->scaley[i]);
803 
804  offsetx_[i] = entry00_->offsetx[i]
805  +adcota_*(entry01_->offsetx[i] - entry00_->offsetx[i])
806  +adcotb_*(entry10_->offsetx[i] - entry00_->offsetx[i]);
807  if(flip_x_) offsetx_[i] = -offsetx_[i];
808 
809  offsety_[i] = entry00_->offsety[i]
810  +adcota_*(entry01_->offsety[i] - entry00_->offsety[i])
811  +adcotb_*(entry10_->offsety[i] - entry00_->offsety[i]);
812  if(flip_y_) offsety_[i] = -offsety_[i];
813  }
814 
815  for(int i=0; i<2 ; ++i) {
816  for(int j=0; j<5 ; ++j) {
817  // Charge loss switches sides when cot(beta) changes sign
818  if(flip_y_) {
819  xypary0x0_[1-i][j] = (float)entry00_->xypar[i][j];
820  xypary1x0_[1-i][j] = (float)entry10_->xypar[i][j];
821  xypary0x1_[1-i][j] = (float)entry01_->xypar[i][j];
822  lanpar_[1-i][j] = entry00_->lanpar[i][j]
823  +adcota_*(entry01_->lanpar[i][j] - entry00_->lanpar[i][j])
824  +adcotb_*(entry10_->lanpar[i][j] - entry00_->lanpar[i][j]);
825  } else {
826  xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
827  xypary1x0_[i][j] = (float)entry10_->xypar[i][j];
828  xypary0x1_[i][j] = (float)entry01_->xypar[i][j];
829  lanpar_[i][j] = entry00_->lanpar[i][j]
830  +adcota_*(entry01_->lanpar[i][j] - entry00_->lanpar[i][j])
831  +adcotb_*(entry10_->lanpar[i][j] - entry00_->lanpar[i][j]);
832  }
833  }
834  }
835 
836  return success_;
837 } // interpolate
838 
839 
840 // *************************************************************************************************************************************
846 // *************************************************************************************************************************************
847 
848 void SiPixelTemplate2D::sideload(SiPixelTemplateEntry2D* entry, int iDtype, float locBx, float locBz, float lorwdy, float lorwdx, float q50, float fbin[3], float xsize, float ysize, float zsize)
849 {
850  // Set class variables to the input parameters
851 
852  entry00_ = entry;
853  entry01_ = entry;
854  entry10_ = entry;
855  Dtype_ = iDtype;
856  lorywidth_ = lorwdy;
857  lorxwidth_ = lorwdx;
858  xsize_ = xsize;
859  ysize_ = ysize;
860  zsize_ = zsize;
861  s50_ = q50;
862  qscale_ = 1.f;
863  for(int i=0; i<3; ++i) {fbin_[i] = fbin[i];}
864 
865  // Set other class variables
866 
867  adcota_ = 0.f;
868  adcotb_ = 0.f;
869 
870  // Interpolate things in cot(alpha)-cot(beta)
871 
872  qavg_ = entry00_->qavg;
873 
874  pixmax_ = entry00_->pixmax;
875 
876  sxymax_ = entry00_->sxymax;
877 
878  clsleny_ = entry00_->clsleny;
879 
880  clslenx_ = entry00_->clslenx;
881 
882  scaleyavg_ = 1.f;
883 
884  scalexavg_ = 1.f;
885 
886  delyavg_ = 0.f;
887 
888  delysig_ = 0.f;
889 
890  for(int i=0; i<4 ; ++i) {
891  scalex_[i] = 1.f;
892  scaley_[i] = 1.f;
893  offsetx_[i] = 0.f;
894  offsety_[i] = 0.f;
895  }
896 
897  // This works only for IP-related tracks
898 
899  flip_x_ = false;
900  flip_y_ = false;
901  float cotbeta = entry00_->cotbeta;
902  switch(Dtype_) {
903  case 0:
904  if(cotbeta < 0.f) {flip_y_ = true;}
905  break;
906  case 1:
907  if(locBz > 0.f) {
908  flip_y_ = true;
909  }
910  break;
911  case 2:
912  case 3:
913  case 4:
914  case 5:
915  if(locBx*locBz < 0.f) {
916  flip_x_ = true;
917  }
918  if(locBx < 0.f) {
919  flip_y_ = true;
920  }
921  break;
922  default:
923 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
924  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::illegal subdetector ID = " << iDtype << std::endl;
925 #else
926  std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << iDtype << std::endl;
927 #endif
928  }
929 
930  // Calculate signed quantities
931 
932  lorydrift_ = lorywidth_/2.;
933  if(flip_y_) lorydrift_ = -lorydrift_;
934  lorxdrift_ = lorxwidth_/2.;
935  if(flip_x_) lorxdrift_ = -lorxdrift_;
936 
937  for(int i=0; i<2 ; ++i) {
938  for(int j=0; j<5 ; ++j) {
939  // Charge loss switches sides when cot(beta) changes sign
940  if(flip_y_) {
941  xypary0x0_[1-i][j] = (float)entry00_->xypar[i][j];
942  xypary1x0_[1-i][j] = (float)entry00_->xypar[i][j];
943  xypary0x1_[1-i][j] = (float)entry00_->xypar[i][j];
944  lanpar_[1-i][j] = entry00_->lanpar[i][j];
945  } else {
946  xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
947  xypary1x0_[i][j] = (float)entry00_->xypar[i][j];
948  xypary0x1_[i][j] = (float)entry00_->xypar[i][j];
949  lanpar_[i][j] = entry00_->lanpar[i][j];
950  }
951  }
952  }
953  return;
954 }
955 
956 
957 // *************************************************************************************************************************************
963 // *************************************************************************************************************************************
964 
965 bool SiPixelTemplate2D::xytemp(float xhit, float yhit, bool ydouble[BYM2], bool xdouble[BXM2], float template2d[BXM2][BYM2], bool derivatives, float dpdx2d[2][BXM2][BYM2], float& QTemplate)
966 {
967  // Interpolate for a new set of track angles
968 
969  // Local variables
970  int pixx, pixy, k0, k1, l0, l1, deltax, deltay, iflipy, jflipx, imin, imax, jmin, jmax;
971  int m, n;
972  float dx, dy, ddx, ddy, adx, ady;
973  // const float deltaxy[2] = {8.33f, 12.5f};
974  const float deltaxy[2] = {16.67f, 25.0f};
975 
976  // Check to see if interpolation is valid
977 
978 
979  // next, determine the indices of the closest point in k (y-displacement), l (x-displacement)
980  // pixy and pixx are the indices of the struck pixel in the (Ty,Tx) system
981  // k0,k1 are the k-indices of the closest and next closest point
982  // l0,l1 are the l-indices of the closest and next closest point
983 
984  pixy = (int)floorf(yhit/ysize_);
985  dy = yhit-(pixy+0.5f)*ysize_;
986  if(flip_y_) {dy = -dy;}
987  k0 = (int)(dy/ysize_*6.f+3.5f);
988  if(k0 < 0) k0 = 0;
989  if(k0 > 6) k0 = 6;
990  ddy = 6.f*dy/ysize_ - (k0-3);
991  ady = fabs(ddy);
992  if(ddy > 0.f) {k1 = k0 + 1; if(k1 > 6) k1 = k0-1;} else {k1 = k0 - 1; if(k1 < 0) k1 = k0+1;}
993  pixx = (int)floorf(xhit/xsize_);
994  dx = xhit-(pixx+0.5f)*xsize_;
995  if(flip_x_) {dx = -dx;}
996  l0 = (int)(dx/xsize_*6.f+3.5f);
997  if(l0 < 0) l0 = 0;
998  if(l0 > 6) l0 = 6;
999  ddx = 6.f*dx/xsize_ - (l0-3);
1000  adx = fabs(ddx);
1001  if(ddx > 0.f) {l1 = l0 + 1; if(l1 > 6) l1 = l0-1;} else {l1 = l0 - 1; if(l1 < 0) l1 = l0+1;}
1002 
1003  // OK, lets do the template interpolation.
1004 
1005  // First find the limits of the indices for non-zero pixels
1006 
1007  imin = std::min(entry00_->iymin,entry10_->iymin);
1008  imin_ = std::min(imin,entry01_->iymin);
1009 
1010  jmin = std::min(entry00_->jxmin,entry10_->jxmin);
1011  jmin_ = std::min(jmin,entry01_->jxmin);
1012 
1013  imax = std::max(entry00_->iymax,entry10_->iymax);
1014  imax_ = std::max(imax,entry01_->iymax);
1015 
1016  jmax = std::max(entry00_->jxmax,entry10_->jxmax);
1017  jmax_ = std::max(jmax,entry01_->jxmax);
1018 
1019  // Calculate the x and y offsets to make the new template
1020 
1021  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1022 
1023  ++pixy; ++pixx;
1024 
1025  // In the template store, the struck pixel is always (THy,THx)
1026 
1027  deltax = pixx - T2HX;
1028  deltay = pixy - T2HY;
1029 
1030  // First zero the local 2-d template
1031 
1032  for(int j=0; j<BXM2; ++j) {for(int i=0; i<BYM2; ++i) {xytemp_[j][i] = 0.f;}}
1033 
1034  // Loop over the non-zero part of the template index space and interpolate
1035 
1036  for(int j=jmin_; j<=jmax_; ++j) {
1037  // Flip indices as needed
1038  if(flip_x_) {jflipx=T2XSIZE-1-j; m = deltax+jflipx;} else {m = deltax+j;}
1039  for(int i=imin_; i<=imax_; ++i) {
1040  if(flip_y_) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
1041  if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
1042  xytemp_[m][n] = (float)entry00_->xytemp[k0][l0][i][j]
1043  + adx*(float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j])
1044  + ady*(float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1045  + adcota_*(float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1046  + adcotb_*(float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1047  }
1048  }
1049  }
1050 
1051  //combine rows and columns to simulate double pixels
1052 
1053  for(int n=1; n<BYM3; ++n) {
1054  if(ydouble[n]) {
1055  // Combine the y-columns
1056  for(int m=1; m<BXM3; ++m) {
1057  xytemp_[m][n] += xytemp_[m][n+1];
1058  }
1059  // Now shift the remaining pixels over by one column
1060  for(int i=n+1; i<BYM3; ++i) {
1061  for(int m=1; m<BXM3; ++m) {
1062  xytemp_[m][i] = xytemp_[m][i+1];
1063  }
1064  }
1065  }
1066  }
1067 
1068  //combine rows and columns to simulate double pixels
1069 
1070  for(int m=1; m<BXM3; ++m) {
1071  if(xdouble[m]) {
1072  // Combine the x-rows
1073  for(int n=1; n<BYM3; ++n) {
1074  xytemp_[m][n] += xytemp_[m+1][n];
1075  }
1076  // Now shift the remaining pixels over by one row
1077  for(int j=m+1; j<BXM3; ++j) {
1078  for(n=1; n<BYM3; ++n) {
1079  xytemp_[j][n] = xytemp_[j+1][n];
1080  }
1081  }
1082  }
1083  }
1084 
1085  // Finally, loop through and increment the external template
1086 
1087  float qtemptot = 0.f;
1088 
1089  for(int n=1; n<BYM3; ++n) {
1090  for(int m=1; m<BXM3; ++m) {
1091  if(xytemp_[m][n] != 0.f) {template2d[m][n] += xytemp_[m][n]; qtemptot += xytemp_[m][n];}
1092  }
1093  }
1094 
1095  QTemplate = qtemptot;
1096 
1097  if(derivatives) {
1098 
1099  float dxytempdx[2][BXM2][BYM2], dxytempdy[2][BXM2][BYM2];
1100 
1101  for(int k = 0; k<2; ++k) {
1102  for(int i = 0; i<BXM2; ++i) {
1103  for(int j = 0; j<BYM2; ++j) {
1104  dxytempdx[k][i][j] = 0.f;
1105  dxytempdy[k][i][j] = 0.f;
1106  dpdx2d[k][i][j] = 0.f;
1107  }
1108  }
1109  }
1110 
1111  // First do shifted +x template
1112 
1113  pixx = (int)floorf((xhit+deltaxy[0])/xsize_);
1114  dx = (xhit+deltaxy[0])-(pixx+0.5f)*xsize_;
1115  if(flip_x_) {dx = -dx;}
1116  l0 = (int)(dx/xsize_*6.f+3.5f);
1117  if(l0 < 0) l0 = 0;
1118  if(l0 > 6) l0 = 6;
1119  ddx = 6.f*dx/xsize_ - (l0-3);
1120  adx = fabs(ddx);
1121  if(ddx > 0.f) {l1 = l0 + 1; if(l1 > 6) l1 = l0-1;} else {l1 = l0 - 1; if(l1 < 0) l1 = l0+1;}
1122 
1123  // OK, lets do the template interpolation.
1124 
1125  // Calculate the x and y offsets to make the new template
1126 
1127  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1128 
1129  ++pixx;
1130 
1131  // In the template store, the struck pixel is always (THy,THx)
1132 
1133  deltax = pixx - T2HX;
1134 
1135  // Loop over the non-zero part of the template index space and interpolate
1136 
1137  for(int j=jmin_; j<=jmax_; ++j) {
1138  // Flip indices as needed
1139  if(flip_x_) {jflipx=T2XSIZE-1-j; m = deltax+jflipx;} else {m = deltax+j;}
1140  for(int i=imin_; i<=imax_; ++i) {
1141  if(flip_y_) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
1142  if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
1143  dxytempdx[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j]
1144  + adx*(float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j])
1145  + ady*(float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1146  + adcota_*(float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1147  + adcotb_*(float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1148  }
1149  }
1150  }
1151 
1152  //combine rows and columns to simulate double pixels
1153 
1154  for(int n=1; n<BYM3; ++n) {
1155  if(ydouble[n]) {
1156  // Combine the y-columns
1157  for(int m=1; m<BXM3; ++m) {
1158  dxytempdx[1][m][n] += dxytempdx[1][m][n+1];
1159  }
1160  // Now shift the remaining pixels over by one column
1161  for(int i=n+1; i<BYM3; ++i) {
1162  for(int m=1; m<BXM3; ++m) {
1163  dxytempdx[1][m][i] = dxytempdx[1][m][i+1];
1164  }
1165  }
1166  }
1167  }
1168 
1169  //combine rows and columns to simulate double pixels
1170 
1171  for(int m=1; m<BXM3; ++m) {
1172  if(xdouble[m]) {
1173  // Combine the x-rows
1174  for(int n=1; n<BYM3; ++n) {
1175  dxytempdx[1][m][n] += dxytempdx[1][m+1][n];
1176  }
1177  // Now shift the remaining pixels over by one row
1178  for(int j=m+1; j<BXM3; ++j) {
1179  for(int n=1; n<BYM3; ++n) {
1180  dxytempdx[1][j][n] = dxytempdx[1][j+1][n];
1181  }
1182  }
1183  }
1184  }
1185 
1186  // Next do shifted -x template
1187 
1188  pixx = (int)floorf((xhit-deltaxy[0])/xsize_);
1189  dx = (xhit-deltaxy[0])-(pixx+0.5f)*xsize_;
1190  if(flip_x_) {dx = -dx;}
1191  l0 = (int)(dx/xsize_*6.f+3.5f);
1192  if(l0 < 0) l0 = 0;
1193  if(l0 > 6) l0 = 6;
1194  ddx = 6.f*dx/xsize_ - (l0-3);
1195  adx = fabs(ddx);
1196  if(ddx > 0.f) {l1 = l0 + 1; if(l1 > 6) l1 = l0-1;} else {l1 = l0 - 1; if(l1 < 0) l1 = l0+1;}
1197 
1198  // OK, lets do the template interpolation.
1199 
1200  // Calculate the x and y offsets to make the new template
1201 
1202  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1203 
1204  ++pixx;
1205 
1206  // In the template store, the struck pixel is always (THy,THx)
1207 
1208  deltax = pixx - T2HX;
1209 
1210  // Loop over the non-zero part of the template index space and interpolate
1211 
1212  for(int j=jmin_; j<=jmax_; ++j) {
1213  // Flip indices as needed
1214  if(flip_x_) {jflipx=T2XSIZE-1-j; m = deltax+jflipx;} else {m = deltax+j;}
1215  for(int i=imin_; i<=imax_; ++i) {
1216  if(flip_y_) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
1217  if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
1218  dxytempdx[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j]
1219  + adx*(float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j])
1220  + ady*(float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1221  + adcota_*(float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1222  + adcotb_*(float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1223  }
1224  }
1225  }
1226 
1227  //combine rows and columns to simulate double pixels
1228 
1229  for(int n=1; n<BYM3; ++n) {
1230  if(ydouble[n]) {
1231  // Combine the y-columns
1232  for(int m=1; m<BXM3; ++m) {
1233  dxytempdx[0][m][n] += dxytempdx[0][m][n+1];
1234  }
1235  // Now shift the remaining pixels over by one column
1236  for(int i=n+1; i<BYM3; ++i) {
1237  for(int m=1; m<BXM3; ++m) {
1238  dxytempdx[0][m][i] = dxytempdx[0][m][i+1];
1239  }
1240  }
1241  }
1242  }
1243 
1244  //combine rows and columns to simulate double pixels
1245 
1246  for(int m=1; m<BXM3; ++m) {
1247  if(xdouble[m]) {
1248  // Combine the x-rows
1249  for(int n=1; n<BYM3; ++n) {
1250  dxytempdx[0][m][n] += dxytempdx[0][m+1][n];
1251  }
1252  // Now shift the remaining pixels over by one row
1253  for(int j=m+1; j<BXM3; ++j) {
1254  for(int n=1; n<BYM3; ++n) {
1255  dxytempdx[0][j][n] = dxytempdx[0][j+1][n];
1256  }
1257  }
1258  }
1259  }
1260 
1261 
1262  // Finally, normalize the derivatives and copy the results to the output array
1263 
1264  for(int n=1; n<BYM3; ++n) {
1265  for(int m=1; m<BXM3; ++m) {
1266  dpdx2d[0][m][n] = (dxytempdx[1][m][n] - dxytempdx[0][m][n])/(2.*deltaxy[0]);
1267  }
1268  }
1269 
1270  // Next, do shifted y template
1271 
1272  pixy = (int)floorf((yhit+deltaxy[1])/ysize_);
1273  dy = (yhit+deltaxy[1])-(pixy+0.5f)*ysize_;
1274  if(flip_y_) {dy = -dy;}
1275  k0 = (int)(dy/ysize_*6.f+3.5f);
1276  if(k0 < 0) k0 = 0;
1277  if(k0 > 6) k0 = 6;
1278  ddy = 6.f*dy/ysize_ - (k0-3);
1279  ady = fabs(ddy);
1280  if(ddy > 0.f) {k1 = k0 + 1; if(k1 > 6) k1 = k0-1;} else {k1 = k0 - 1; if(k1 < 0) k1 = k0+1;}
1281  pixx = (int)floorf(xhit/xsize_);
1282  dx = xhit-(pixx+0.5f)*xsize_;
1283  if(flip_x_) {dx = -dx;}
1284  l0 = (int)(dx/xsize_*6.f+3.5f);
1285  if(l0 < 0) l0 = 0;
1286  if(l0 > 6) l0 = 6;
1287  ddx = 6.f*dx/xsize_ - (l0-3);
1288  adx = fabs(ddx);
1289  if(ddx > 0.f) {l1 = l0 + 1; if(l1 > 6) l1 = l0-1;} else {l1 = l0 - 1; if(l1 < 0) l1 = l0+1;}
1290 
1291  // OK, lets do the template interpolation.
1292 
1293  // Calculate the x and y offsets to make the new template
1294 
1295  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1296 
1297  ++pixy; ++pixx;
1298 
1299  // In the template store, the struck pixel is always (THy,THx)
1300 
1301  deltax = pixx - T2HX;
1302  deltay = pixy - T2HY;
1303 
1304  // Loop over the non-zero part of the template index space and interpolate
1305 
1306  for(int j=jmin_; j<=jmax_; ++j) {
1307  // Flip indices as needed
1308  if(flip_x_) {jflipx=T2XSIZE-1-j; m = deltax+jflipx;} else {m = deltax+j;}
1309  for(int i=imin_; i<=imax_; ++i) {
1310  if(flip_y_) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
1311  if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
1312  dxytempdy[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j]
1313  + adx*(float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j])
1314  + ady*(float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1315  + adcota_*(float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1316  + adcotb_*(float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1317  }
1318  }
1319  }
1320 
1321  //combine rows and columns to simulate double pixels
1322 
1323  for(int n=1; n<BYM3; ++n) {
1324  if(ydouble[n]) {
1325  // Combine the y-columns
1326  for(int m=1; m<BXM3; ++m) {
1327  dxytempdy[1][m][n] += dxytempdy[1][m][n+1];
1328  }
1329  // Now shift the remaining pixels over by one column
1330  for(int i=n+1; i<BYM3; ++i) {
1331  for(int m=1; m<BXM3; ++m) {
1332  dxytempdy[1][m][i] = dxytempdy[1][m][i+1];
1333  }
1334  }
1335  }
1336  }
1337 
1338  //combine rows and columns to simulate double pixels
1339 
1340  for(int m=1; m<BXM3; ++m) {
1341  if(xdouble[m]) {
1342  // Combine the x-rows
1343  for(int n=1; n<BYM3; ++n) {
1344  dxytempdy[1][m][n] += dxytempdy[1][m+1][n];
1345  }
1346  // Now shift the remaining pixels over by one row
1347  for(int j=m+1; j<BXM3; ++j) {
1348  for(int n=1; n<BYM3; ++n) {
1349  dxytempdy[1][j][n] = dxytempdy[1][j+1][n];
1350  }
1351  }
1352  }
1353  }
1354 
1355  // Next, do shifted -y template
1356 
1357  pixy = (int)floorf((yhit-deltaxy[1])/ysize_);
1358  dy = (yhit-deltaxy[1])-(pixy+0.5f)*ysize_;
1359  if(flip_y_) {dy = -dy;}
1360  k0 = (int)(dy/ysize_*6.f+3.5f);
1361  if(k0 < 0) k0 = 0;
1362  if(k0 > 6) k0 = 6;
1363  ddy = 6.f*dy/ysize_ - (k0-3);
1364  ady = fabs(ddy);
1365  if(ddy > 0.f) {k1 = k0 + 1; if(k1 > 6) k1 = k0-1;} else {k1 = k0 - 1; if(k1 < 0) k1 = k0+1;}
1366 
1367  // OK, lets do the template interpolation.
1368 
1369  // Calculate the x and y offsets to make the new template
1370 
1371  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1372 
1373  ++pixy;
1374 
1375  // In the template store, the struck pixel is always (THy,THx)
1376 
1377  deltay = pixy - T2HY;
1378 
1379  // Loop over the non-zero part of the template index space and interpolate
1380 
1381  for(int j=jmin_; j<=jmax_; ++j) {
1382  // Flip indices as needed
1383  if(flip_x_) {jflipx=T2XSIZE-1-j; m = deltax+jflipx;} else {m = deltax+j;}
1384  for(int i=imin_; i<=imax_; ++i) {
1385  if(flip_y_) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
1386  if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
1387  dxytempdy[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j]
1388  + adx*(float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j])
1389  + ady*(float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1390  + adcota_*(float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1391  + adcotb_*(float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1392  }
1393  }
1394  }
1395 
1396  //combine rows and columns to simulate double pixels
1397 
1398  for(int n=1; n<BYM3; ++n) {
1399  if(ydouble[n]) {
1400  // Combine the y-columns
1401  for(int m=1; m<BXM3; ++m) {
1402  dxytempdy[0][m][n] += dxytempdy[0][m][n+1];
1403  }
1404  // Now shift the remaining pixels over by one column
1405  for(int i=n+1; i<BYM3; ++i) {
1406  for(int m=1; m<BXM3; ++m) {
1407  dxytempdy[0][m][i] = dxytempdy[0][m][i+1];
1408  }
1409  }
1410  }
1411  }
1412 
1413  //combine rows and columns to simulate double pixels
1414 
1415  for(int m=1; m<BXM3; ++m) {
1416  if(xdouble[m]) {
1417  // Combine the x-rows
1418  for(int n=1; n<BYM3; ++n) {
1419  dxytempdy[0][m][n] += dxytempdy[0][m+1][n];
1420  }
1421  // Now shift the remaining pixels over by one row
1422  for(int j=m+1; j<BXM3; ++j) {
1423  for(int n=1; n<BYM3; ++n) {
1424  dxytempdy[0][j][n] = dxytempdy[0][j+1][n];
1425  }
1426  }
1427  }
1428  }
1429 
1430  // Finally, normalize the derivatives and copy the results to the output array
1431 
1432  for(int n=1; n<BYM3; ++n) {
1433  for(int m=1; m<BXM3; ++m) {
1434  dpdx2d[1][m][n] = (dxytempdy[1][m][n] - dxytempdy[0][m][n])/(2.*deltaxy[1]);
1435  }
1436  }
1437  }
1438 
1439  return success_;
1440 } // xytemp
1441 
1442 
1443 
1444 // *************************************************************************************************************************************
1451 // *************************************************************************************************************************************
1452 
1453 bool SiPixelTemplate2D::xytemp(float xhit, float yhit, bool ydouble[BYM2], bool xdouble[BXM2], float template2d[BXM2][BYM2])
1454 {
1455  // Interpolate for a new set of track angles
1456 
1457  bool derivatives = false;
1458  float dpdx2d[2][BXM2][BYM2];
1459  float QTemplate;
1460 
1461  return SiPixelTemplate2D::xytemp(xhit, yhit, ydouble, xdouble, template2d, derivatives, dpdx2d, QTemplate);
1462 
1463 } // xytemp
1464 
1465 
1466 
1467 // *************************************************************************************************************************************
1477 // *************************************************************************************************************************************
1478 
1479 bool SiPixelTemplate2D::xytemp(int id, float cotalpha, float cotbeta, float xhit, float yhit, std::vector<bool>& ydouble, std::vector<bool>& xdouble, float template2d[BXM2][BYM2])
1480 {
1481 
1482  // Local variables
1483 
1484  bool derivatives = false;
1485  float dpdx2d[2][BXM2][BYM2];
1486  float QTemplate;
1487  float locBx = 1.f;
1488  if(cotbeta < 0.f) {locBx = -1.f;}
1489  float locBz = locBx;
1490  if(cotalpha < 0.f) {locBz = -locBx;}
1491 
1492  bool yd[BYM2], xd[BXM2];
1493 
1494  yd[0] = false; yd[BYM2-1] = false;
1495  for(int i=0; i<TYSIZE; ++i) { yd[i+1] = ydouble[i];}
1496  xd[0] = false; xd[BXM2-1] = false;
1497  for(int j=0; j<TXSIZE; ++j) { xd[j+1] = xdouble[j];}
1498 
1499 
1500  // Interpolate for a new set of track angles
1501 
1502  if(SiPixelTemplate2D::interpolate(id, cotalpha, cotbeta, locBz, locBx)) {
1503  return SiPixelTemplate2D::xytemp(xhit, yhit, yd, xd, template2d, derivatives, dpdx2d, QTemplate);
1504  } else {
1505  return false;
1506  }
1507 
1508 } // xytemp
1509 
1510 
1511 // ************************************************************************************************************
1517 // ************************************************************************************************************
1518 void SiPixelTemplate2D::xysigma2(float qpixel, int index, float& xysig2)
1519 
1520 {
1521  // Interpolate using quantities already stored in the private variables
1522 
1523  // Local variables
1524  float sigi, sigi2, sigi3, sigi4, qscale, err2, err00;
1525 
1526  // Make sure that input is OK
1527 
1528 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1529  if(index < 1 || index >= BYM2) {
1530  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::ysigma2 called with index = " << index << std::endl;
1531  }
1532 #else
1533  assert(index > 0 && index < BYM2);
1534 #endif
1535 
1536  // Define the maximum signal to use in the parameterization
1537 
1538  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1539 
1540  if(qpixel < sxymax_) {
1541  sigi = qpixel;
1542  qscale = 1.f;
1543  } else {
1544  sigi = sxymax_;
1545  qscale = qpixel/sxymax_;
1546  }
1547  sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
1548  if(index <= T2HYP1) {
1549  err00 = xypary0x0_[0][0]+xypary0x0_[0][1]*sigi+xypary0x0_[0][2]*sigi2+xypary0x0_[0][3]*sigi3+xypary0x0_[0][4]*sigi4;
1550  err2 = err00
1551  +adcota_*(xypary0x1_[0][0]+xypary0x1_[0][1]*sigi+xypary0x1_[0][2]*sigi2+xypary0x1_[0][3]*sigi3+xypary0x1_[0][4]*sigi4 - err00)
1552  +adcotb_*(xypary1x0_[0][0]+xypary1x0_[0][1]*sigi+xypary1x0_[0][2]*sigi2+xypary1x0_[0][3]*sigi3+xypary1x0_[0][4]*sigi4 - err00);
1553  } else {
1554  err00 = xypary0x0_[1][0]+xypary0x0_[1][1]*sigi+xypary0x0_[1][2]*sigi2+xypary0x0_[1][3]*sigi3+xypary0x0_[1][4]*sigi4;
1555  err2 = err00
1556  +adcota_*(xypary0x1_[1][0]+xypary0x1_[1][1]*sigi+xypary0x1_[1][2]*sigi2+xypary0x1_[1][3]*sigi3+xypary0x1_[1][4]*sigi4 - err00)
1557  +adcotb_*(xypary1x0_[1][0]+xypary1x0_[1][1]*sigi+xypary1x0_[1][2]*sigi2+xypary1x0_[1][3]*sigi3+xypary1x0_[1][4]*sigi4 - err00);
1558  }
1559  xysig2 =qscale*err2;
1560  if(xysig2 <= 0.f) {xysig2 = s50_*s50_;}
1561 
1562  return;
1563 
1564 } // End xysigma2
1565 
1566 
1567 // ************************************************************************************************************
1569 // ************************************************************************************************************
1570 void SiPixelTemplate2D::landau_par(float lanpar[2][5])
1571 
1572 {
1573  // Interpolate using quantities already stored in the private variables
1574 
1575  for(int i=0; i<2; ++i) {
1576  for(int j=0; j<5; ++j) {
1577  lanpar[i][j] = lanpar_[i][j];
1578  }
1579  }
1580  return;
1581 
1582 } // End lan_par
1583 
1584 
1585 
SiPixelTemplateEntry2D::pixmax
float pixmax
maximum charge for individual pixels in cluster
Definition: SiPixelTemplate2D.h:48
SiPixelTemplateHeader2D::qscale
float qscale
Charge scaling to match cmssw and pixelav.
Definition: SiPixelTemplate2D.h:86
T2HYP1
#define T2HYP1
Definition: SiPixelTemplateDefs.h:41
mps_fire.i
i
Definition: mps_fire.py:338
SiPixelTemplateEntry2D::qavg
float qavg
average cluster charge for this set of track angles
Definition: SiPixelTemplate2D.h:47
SiPixelTemplateHeader2D::Vbias
float Vbias
detector bias potential in Volts
Definition: SiPixelTemplate2D.h:91
SiPixelTemplateHeader2D::lorxwidth
float lorxwidth
estimate of x-lorentz width for optimal resolution
Definition: SiPixelTemplate2D.h:88
SiPixelTemplate2D::xytemp
bool xytemp(float xhit, float yhit, bool ydouble[21+2], bool xdouble[13+2], float template2d[13+2][21+2], bool dervatives, float dpdx2d[2][13+2][21+2], float &QTemplate)
Definition: SiPixelTemplate2D.cc:965
SiPixelTemplateEntry2D::runnum
int runnum
< Basic template entry corresponding to a single set of track angles
Definition: SiPixelTemplate2D.h:43
SiPixelTemplateHeader2D::ysize
float ysize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate2D.h:101
SiPixelTemplateStore2D::head
SiPixelTemplateHeader2D head
< template storage structure
Definition: SiPixelTemplate2D.h:108
SiPixelTemplate2D::xysigma2
void xysigma2(float qpixel, int index, float &xysig2)
Definition: SiPixelTemplate2D.cc:1518
TXSIZE
#define TXSIZE
Definition: SiPixelTemplateDefs.h:30
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
T2HY
#define T2HY
Definition: SiPixelTemplateDefs.h:40
SiPixelTemplateEntry2D::clsleny
float clsleny
cluster y-length in pixels at signal height symax/2
Definition: SiPixelTemplate2D.h:61
SiPixelTemplateEntry2D::sxymax
float sxymax
average pixel signal for use of the error parameterization
Definition: SiPixelTemplate2D.h:49
BXM3
#define BXM3
Definition: SiPixelTemplateDefs.h:37
T2HX
#define T2HX
Definition: SiPixelTemplateDefs.h:42
SiPixelTemplateEntry2D::xypar
float xypar[2][5]
pixel uncertainty parameterization
Definition: SiPixelTemplate2D.h:54
SiPixel2DTemplateDBObject::sVector
std::vector< float > sVector() const
Definition: SiPixel2DTemplateDBObject.h:67
ENDL
#define ENDL
Definition: SiPixelTemplate2D.cc:47
SiPixelTemplateEntry2D::scaleyavg
float scaleyavg
average y-error scale factor
Definition: SiPixelTemplate2D.h:67
ysize
const Int_t ysize
Definition: trackSplitPlot.h:43
SiPixelTemplateEntry2D::delyavg
float delyavg
average length difference between template and cluster
Definition: SiPixelTemplate2D.h:68
SiPixelTemplateEntry2D::delysig
float delysig
rms of length difference between template and cluster
Definition: SiPixelTemplate2D.h:69
edm::isFinite
constexpr bool isFinite(T x)
SiPixelTemplateStore2D::entry
SiPixelTemplateEntry2D ** entry
use 2d entry to store BPix and FPix entries [dynamically allocated]
Definition: SiPixelTemplate2D.h:109
SiPixelTemplateEntry2D::cotbeta
float cotbeta
cot(beta) is proportional to cluster length in y and is basis of interpolation
Definition: SiPixelTemplate2D.h:45
BXM2
#define BXM2
Definition: SiPixelTemplateDefs.h:36
SiPixelTemplateHeader2D::lorxbias
float lorxbias
estimate of x-lorentz bias
Definition: SiPixelTemplate2D.h:90
SiPixelTemplateEntry2D::iymax
int iymax
the maximum nonzero pixel yindex in template (saves time during interpolation)
Definition: SiPixelTemplate2D.h:51
SiPixelTemplateHeader2D::fbin
float fbin[3]
The QBin definitions in Q_clus/Q_avg.
Definition: SiPixelTemplate2D.h:99
T2XSIZE
#define T2XSIZE
Definition: SiPixelTemplateDefs.h:39
mps_splice.entry
entry
Definition: mps_splice.py:68
SiPixelTemplateHeader2D::ss50
float ss50
1/2 of the single hit dcol threshold in electrons
Definition: SiPixelTemplate2D.h:95
SiPixelTemplateEntry2D::iymin
int iymin
the minimum nonzero pixel yindex in template (saves time during interpolation)
Definition: SiPixelTemplate2D.h:50
SiPixelTemplateEntry2D::xytemp
short int xytemp[7][7][21][7]
templates for y-reconstruction (binned over 1 central pixel)
Definition: SiPixelTemplate2D.h:56
SiPixelTemplateHeader2D::ID
int ID
< template header structure
Definition: SiPixelTemplate2D.h:81
LOGERROR
#define LOGERROR(x)
Definition: SiPixelTemplate2D.cc:45
triggerObjects_cff.id
id
Definition: triggerObjects_cff.py:29
SiPixelTemplateHeader2D::NTxx
int NTxx
number of Template x-entries in each slice
Definition: SiPixelTemplate2D.h:84
SiPixelTemplateEntry2D::mpvvav
float mpvvav
most probable charge in Vavilov distribution (not actually for larger kappa)
Definition: SiPixelTemplate2D.h:63
LOGINFO
#define LOGINFO(x)
Definition: SiPixelTemplate2D.cc:46
FrontierConditions_GlobalTag_cff.file
file
Definition: FrontierConditions_GlobalTag_cff.py:13
SiPixelTemplateEntry2D::scalex
float scalex[4]
x-error scale factor in 4 charge bins
Definition: SiPixelTemplate2D.h:70
SiPixelTemplateEntry2D::jxmax
int jxmax
the maximum nonzero pixel xindex in template (saves time during interpolation)
Definition: SiPixelTemplate2D.h:53
SiPixelTemplateEntry2D::lanpar
float lanpar[2][5]
pixel landau distribution parameters
Definition: SiPixelTemplate2D.h:55
SiPixelTemplateHeader2D::xsize
float xsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate2D.h:100
SiPixelTemplateEntry2D::costrk
float costrk[3]
direction cosines of tracks used to generate this entry
Definition: SiPixelTemplate2D.h:46
SiPixelTemplateEntry2D::offsetx
float offsetx[4]
x-offset in 4 charge bins
Definition: SiPixelTemplate2D.h:72
BYM2
#define BYM2
Definition: SiPixelTemplateDefs.h:28
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
SiPixelTemplateEntry2D::jxmin
int jxmin
the minimum nonzero pixel xindex in template (saves time during interpolation)
Definition: SiPixelTemplate2D.h:52
min
T min(T a, T b)
Definition: MathUtil.h:58
BYM3
#define BYM3
Definition: SiPixelTemplateDefs.h:29
SiPixelTemplateEntry2D::sigmavav
float sigmavav
"sigma" scale fctor for Vavilov distribution
Definition: SiPixelTemplate2D.h:64
SiPixelTemplateEntry2D::cotalpha
float cotalpha
cot(alpha) is proportional to cluster length in x and is basis of interpolation
Definition: SiPixelTemplate2D.h:44
SiPixelTemplate2D::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
Definition: SiPixelTemplate2D.cc:614
SiPixelTemplateHeader2D::title
char title[80]
template title
Definition: SiPixelTemplate2D.h:96
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
SiPixelTemplateHeader2D::zsize
float zsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelTemplate2D.h:102
SiPixelTemplateEntry2D::offsety
float offsety[4]
y-offset in 4 charge bins
Definition: SiPixelTemplate2D.h:73
TYSIZE
#define TYSIZE
Definition: SiPixelTemplateDefs.h:21
SiPixelTemplateEntry2D::chi2scale
float chi2scale
scale factor for the chi2 distribution
Definition: SiPixelTemplate2D.h:58
SiPixelTemplateHeader2D::lorywidth
float lorywidth
estimate of y-lorentz width for optimal resolution
Definition: SiPixelTemplate2D.h:87
SiPixelTemplateEntry2D::clslenx
float clslenx
cluster x-length in pixels at signal height sxmax/2
Definition: SiPixelTemplate2D.h:62
SiPixelTemplateHeader2D::NTy
int NTy
number of Template y entries
Definition: SiPixelTemplate2D.h:82
SiPixelTemplate2D::pushfile
static bool pushfile(int filenum, std::vector< SiPixelTemplateStore2D > &pixelTemp, std::string dir="")
Definition: SiPixelTemplate2D.cc:63
SiPixelTemplateEntry2D::scalexavg
float scalexavg
average x-error scale factor
Definition: SiPixelTemplate2D.h:66
SiPixelTemplateHeader2D::Dtype
int Dtype
detector type (0=BPix, 1=FPix)
Definition: SiPixelTemplate2D.h:85
SiPixelTemplateEntry2D::scaley
float scaley[4]
y-error scale factor in 4 charge bins
Definition: SiPixelTemplate2D.h:71
edm
HLT enums.
Definition: AlignableModifier.h:17
SiPixelTemplate2D::landau_par
void landau_par(float lanpar[2][5])
Return the Landau probability parameters for this set of cot(alpha, cot(beta)
Definition: SiPixelTemplate2D.cc:1570
edm::FileInPath::fullPath
std::string fullPath() const
Definition: FileInPath.cc:163
SiPixelTemplateEntry2D::chi2ppix
float chi2ppix
average chi^2 per pixel
Definition: SiPixelTemplate2D.h:57
SiPixelTemplateHeader2D::templ_version
int templ_version
Version number of the template to ensure code compatibility.
Definition: SiPixelTemplate2D.h:97
gather_cfg.cout
cout
Definition: gather_cfg.py:144
dir
dbl *** dir
Definition: mlp_gen.cc:35
SiPixelTemplateHeader2D::s50
float s50
1/2 of the multihit dcol threshold in electrons
Definition: SiPixelTemplate2D.h:94
SiPixelTemplateHeader2D::temperature
float temperature
detector temperature in deg K
Definition: SiPixelTemplate2D.h:92
SiPixelTemplateHeader2D::Bfield
float Bfield
Bfield in Tesla.
Definition: SiPixelTemplate2D.h:98
SiPixelTemplate2D::sideload
void sideload(SiPixelTemplateEntry2D *entry, int iDtype, float locBx, float locBz, float lorwdy, float lorwdx, float q50, float fbin[3], float xsize, float ysize, float zsize)
Definition: SiPixelTemplate2D.cc:848
T2YSIZE
#define T2YSIZE
Definition: SiPixelTemplateDefs.h:38
xsize
const Int_t xsize
Definition: trackSplitPlot.h:42
SiPixelTemplateHeader2D::fluence
float fluence
radiation fluence in n_eq/cm^2
Definition: SiPixelTemplate2D.h:93
SiPixelTemplateEntry2D::kappavav
float kappavav
kappa parameter for Vavilov distribution
Definition: SiPixelTemplate2D.h:65
SiPixelTemplateHeader2D::lorybias
float lorybias
estimate of y-lorentz bias
Definition: SiPixelTemplate2D.h:89
SiPixelTemplateHeader2D::NTyx
int NTyx
number of Template y-slices of x entries
Definition: SiPixelTemplate2D.h:83
Generated for CMSSW Reference Manual by doxygen 1.8.11