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