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 bool SiPixelTemplate2D::getid(int id)
523 {
524  if(id != id_current_) {
525 
526  // Find the index corresponding to id
527 
528  index_id_ = -1;
529  for(int i=0; i<(int)thePixelTemp_.size(); ++i) {
530 
531  if(id == thePixelTemp_[i].head.ID) {
532 
533  index_id_ = i;
534  id_current_ = id;
535 
536  // Copy the charge scaling factor to the private variable
537 
538  Dtype_ = thePixelTemp_[index_id_].head.Dtype;
539 
540  // Copy the charge scaling factor to the private variable
541 
542  qscale_ = thePixelTemp_[index_id_].head.qscale;
543 
544  // Copy the pseudopixel signal size to the private variable
545 
546  s50_ = thePixelTemp_[index_id_].head.s50;
547 
548  // Copy Qbinning info to private variables
549 
550  for(int j=0; j<3; ++j) {fbin_[j] = thePixelTemp_[index_id_].head.fbin[j];}
551 
552  // Copy the Lorentz widths to private variables
553 
554  lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
555  lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
556 
557  // Copy the pixel sizes private variables
558 
559  xsize_ = thePixelTemp_[index_id_].head.xsize;
560  ysize_ = thePixelTemp_[index_id_].head.ysize;
561  zsize_ = thePixelTemp_[index_id_].head.zsize;
562 
563  // Determine the size of this template
564 
565  Nyx_ = thePixelTemp_[index_id_].head.NTyx;
566  Nxx_ = thePixelTemp_[index_id_].head.NTxx;
567 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
568  if(Nyx_ < 2 || Nxx_ < 2) {
569  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Nyx/Nxx = " << Nyx_ << "/" << Nxx_ << std::endl;
570  }
571 #else
572  assert(Nyx_ > 1 && Nxx_ > 1);
573 #endif
574  int imidx = Nxx_/2;
575 
576  cotalpha0_ = thePixelTemp_[index_id_].entry[0][0].cotalpha;
577  cotalpha1_ = thePixelTemp_[index_id_].entry[0][Nxx_-1].cotalpha;
578  deltacota_ = (cotalpha1_-cotalpha0_)/(float)(Nxx_-1);
579 
580  cotbeta0_ = thePixelTemp_[index_id_].entry[0][imidx].cotbeta;
581  cotbeta1_ = thePixelTemp_[index_id_].entry[Nyx_-1][imidx].cotbeta;
582  deltacotb_ = (cotbeta1_-cotbeta0_)/(float)(Nyx_-1);
583 
584  break;
585  }
586  }
587  }
588 
589 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
590  if(index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
591  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
592  << ", Are you using the correct global tag?" << std::endl;
593  }
594 #else
595  assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
596 #endif
597  return true;
598 }
599 
600 // *************************************************************************************************************************************
611 // *************************************************************************************************************************************
612 
613 bool SiPixelTemplate2D::interpolate(int id, float cotalpha, float cotbeta, float locBz, float locBx)
614 {
615 
616 
617  // Interpolate for a new set of track angles
618 
619  // Local variables
620 
621  float acotb, dcota, dcotb;
622 
623  // Check to see if interpolation is valid
624 
625  if(id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
626  cota_current_ = cotalpha; cotb_current_ = cotbeta;
627  success_ = getid(id);
628  }
629 
630 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
631  if(index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
632  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
633  << ", Are you using the correct global tag?" << std::endl;
634  }
635 #else
636  assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
637 #endif
638 
639  // Check angle limits and et up interpolation parameters
640 
641  float cota = cotalpha;
642  flip_x_ = false;
643  flip_y_ = false;
644  switch(Dtype_) {
645  case 0:
646  if(cotbeta < 0.f) {flip_y_ = true;}
647  break;
648  case 1:
649  if(locBz > 0.f) {flip_y_ = true;}
650  break;
651  case 2:
652  case 3:
653  case 4:
654  case 5:
655  if(locBx*locBz < 0.f) {
656  cota = fabs(cotalpha);
657  flip_x_ = true;
658  }
659  if(locBx < 0.f) {flip_y_ = true;}
660  break;
661  default:
662 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
663  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
664 #else
665  std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << thePixelTemp_[index_id_].head.Dtype << std::endl;
666 #endif
667  }
668 
669  if(cota < cotalpha0_) {
670  success_ = false;
671  jx0_ = 0;
672  jx1_ = 1;
673  adcota_ = 0.f;
674  } else if(cota > cotalpha1_) {
675  success_ = false;
676  jx0_ = Nxx_ - 1;
677  jx1_ = jx0_ - 1;
678  adcota_ = 0.f;
679  } else {
680  jx0_ = (int)((cota-cotalpha0_)/deltacota_+0.5f);
681  dcota = (cota - (cotalpha0_ + jx0_*deltacota_))/deltacota_;
682  adcota_ = fabs(dcota);
683  if(dcota > 0.f) {jx1_ = jx0_ + 1;if(jx1_ > Nxx_-1) jx1_ = jx0_-1;} else {jx1_ = jx0_ - 1; if(jx1_ < 0) jx1_ = jx0_+1;}
684  }
685 
686  // Interpolate the absolute value of cot(beta)
687 
688  acotb = fabs(cotbeta);
689 
690  if(acotb < cotbeta0_) {
691  success_ = false;
692  iy0_ = 0;
693  iy1_ = 1;
694  adcotb_ = 0.f;
695  } else if(acotb > cotbeta1_) {
696  success_ = false;
697  iy0_ = Nyx_ - 1;
698  iy1_ = iy0_ - 1;
699  adcotb_ = 0.f;
700  } else {
701  iy0_ = (int)((acotb-cotbeta0_)/deltacotb_+0.5f);
702  dcotb = (acotb - (cotbeta0_ + iy0_*deltacotb_))/deltacotb_;
703  adcotb_ = fabs(dcotb);
704  if(dcotb > 0.f) {iy1_ = iy0_ + 1; if(iy1_ > Nyx_-1) iy1_ = iy0_-1;} else {iy1_ = iy0_ - 1; if(iy1_ < 0) iy1_ = iy0_+1;}
705  }
706 
707  // Calculate signed quantities
708 
709  lorydrift_ = lorywidth_/2.;
710  if(flip_y_) lorydrift_ = -lorydrift_;
711  lorxdrift_ = lorxwidth_/2.;
712  if(flip_x_) lorxdrift_ = -lorxdrift_;
713 
714  // Use pointers to the three angle pairs used in the interpolation
715 
716 
717  entry00_ = &thePixelTemp_[index_id_].entry[iy0_][jx0_];
718  entry10_ = &thePixelTemp_[index_id_].entry[iy1_][jx0_];
719  entry01_ = &thePixelTemp_[index_id_].entry[iy0_][jx1_];
720 
721  // Interpolate things in cot(alpha)-cot(beta)
722 
723  qavg_ = entry00_->qavg
724  +adcota_*(entry01_->qavg - entry00_->qavg)
725  +adcotb_*(entry10_->qavg - entry00_->qavg);
726 
727  pixmax_ = entry00_->pixmax
728  +adcota_*(entry01_->pixmax - entry00_->pixmax)
729  +adcotb_*(entry10_->pixmax - entry00_->pixmax);
730 
731  sxymax_ = entry00_->sxymax
732  +adcota_*(entry01_->sxymax - entry00_->sxymax)
733  +adcotb_*(entry10_->sxymax - entry00_->sxymax);
734 
735  chi2avgone_ = entry00_->chi2avgone
736  +adcota_*(entry01_->chi2avgone - entry00_->chi2avgone)
737  +adcotb_*(entry10_->chi2avgone - entry00_->chi2avgone);
738 
739  chi2minone_ = entry00_->chi2minone
740  +adcota_*(entry01_->chi2minone - entry00_->chi2minone)
741  +adcotb_*(entry10_->chi2minone - entry00_->chi2minone);
742 
743  clsleny_ = entry00_->clsleny
744  +adcota_*(entry01_->clsleny - entry00_->clsleny)
745  +adcotb_*(entry10_->clsleny - entry00_->clsleny);
746 
747  clslenx_ = entry00_->clslenx
748  +adcota_*(entry01_->clslenx - entry00_->clslenx)
749  +adcotb_*(entry10_->clslenx - entry00_->clslenx);
750 
751 
752  chi2ppix_ = entry00_->chi2ppix
753  +adcota_*(entry01_->chi2ppix - entry00_->chi2ppix)
754  +adcotb_*(entry10_->chi2ppix - entry00_->chi2ppix);
755 
756  chi2scale_ = entry00_->chi2scale
757  +adcota_*(entry01_->chi2scale - entry00_->chi2scale)
758  +adcotb_*(entry10_->chi2scale - entry00_->chi2scale);
759 
760  scaleyavg_ = entry00_->scaleyavg
761  +adcota_*(entry01_->scaleyavg - entry00_->scaleyavg)
762  +adcotb_*(entry10_->scaleyavg - entry00_->scaleyavg);
763 
764  scalexavg_ = entry00_->scalexavg
765  +adcota_*(entry01_->scalexavg - entry00_->scalexavg)
766  +adcotb_*(entry10_->scalexavg - entry00_->scalexavg);
767 
768  delyavg_ = entry00_->delyavg
769  +adcota_*(entry01_->delyavg - entry00_->delyavg)
770  +adcotb_*(entry10_->delyavg - entry00_->delyavg);
771 
772  delysig_ = entry00_->delysig
773  +adcota_*(entry01_->delysig - entry00_->delysig)
774  +adcotb_*(entry10_->delysig - entry00_->delysig);
775 
776  mpvvav_ = entry00_->mpvvav
777  +adcota_*(entry01_->mpvvav - entry00_->mpvvav)
778  +adcotb_*(entry10_->mpvvav - entry00_->mpvvav);
779 
780  sigmavav_ = entry00_->sigmavav
781  +adcota_*(entry01_->sigmavav - entry00_->sigmavav)
782  +adcotb_*(entry10_->sigmavav - entry00_->sigmavav);
783 
784  kappavav_ = entry00_->kappavav
785  +adcota_*(entry01_->kappavav - entry00_->kappavav)
786  +adcotb_*(entry10_->kappavav - entry00_->kappavav);
787 
788  for(int i=0; i<4 ; ++i) {
789  scalex_[i] = entry00_->scalex[i]
790  +adcota_*(entry01_->scalex[i] - entry00_->scalex[i])
791  +adcotb_*(entry10_->scalex[i] - entry00_->scalex[i]);
792 
793  scaley_[i] = entry00_->scaley[i]
794  +adcota_*(entry01_->scaley[i] - entry00_->scaley[i])
795  +adcotb_*(entry10_->scaley[i] - entry00_->scaley[i]);
796 
797  offsetx_[i] = entry00_->offsetx[i]
798  +adcota_*(entry01_->offsetx[i] - entry00_->offsetx[i])
799  +adcotb_*(entry10_->offsetx[i] - entry00_->offsetx[i]);
800  if(flip_x_) offsetx_[i] = -offsetx_[i];
801 
802  offsety_[i] = entry00_->offsety[i]
803  +adcota_*(entry01_->offsety[i] - entry00_->offsety[i])
804  +adcotb_*(entry10_->offsety[i] - entry00_->offsety[i]);
805  if(flip_y_) offsety_[i] = -offsety_[i];
806  }
807 
808  for(int i=0; i<2 ; ++i) {
809  for(int j=0; j<5 ; ++j) {
810  // Charge loss switches sides when cot(beta) changes sign
811  if(flip_y_) {
812  xypary0x0_[1-i][j] = (float)entry00_->xypar[i][j];
813  xypary1x0_[1-i][j] = (float)entry10_->xypar[i][j];
814  xypary0x1_[1-i][j] = (float)entry01_->xypar[i][j];
815  lanpar_[1-i][j] = entry00_->lanpar[i][j]
816  +adcota_*(entry01_->lanpar[i][j] - entry00_->lanpar[i][j])
817  +adcotb_*(entry10_->lanpar[i][j] - entry00_->lanpar[i][j]);
818  } else {
819  xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
820  xypary1x0_[i][j] = (float)entry10_->xypar[i][j];
821  xypary0x1_[i][j] = (float)entry01_->xypar[i][j];
822  lanpar_[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  }
826  }
827  }
828 
829  return success_;
830 } // interpolate
831 
832 
833 // *************************************************************************************************************************************
839 // *************************************************************************************************************************************
840 
841 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)
842 {
843  // Set class variables to the input parameters
844 
845  entry00_ = entry;
846  entry01_ = entry;
847  entry10_ = entry;
848  Dtype_ = iDtype;
849  lorywidth_ = lorwdy;
850  lorxwidth_ = lorwdx;
851  xsize_ = xsize;
852  ysize_ = ysize;
853  zsize_ = zsize;
854  s50_ = q50;
855  qscale_ = 1.f;
856  for(int i=0; i<3; ++i) {fbin_[i] = fbin[i];}
857 
858  // Set other class variables
859 
860  adcota_ = 0.f;
861  adcotb_ = 0.f;
862 
863  // Interpolate things in cot(alpha)-cot(beta)
864 
865  qavg_ = entry00_->qavg;
866 
867  pixmax_ = entry00_->pixmax;
868 
869  sxymax_ = entry00_->sxymax;
870 
871  clsleny_ = entry00_->clsleny;
872 
873  clslenx_ = entry00_->clslenx;
874 
875  scaleyavg_ = 1.f;
876 
877  scalexavg_ = 1.f;
878 
879  delyavg_ = 0.f;
880 
881  delysig_ = 0.f;
882 
883  for(int i=0; i<4 ; ++i) {
884  scalex_[i] = 1.f;
885  scaley_[i] = 1.f;
886  offsetx_[i] = 0.f;
887  offsety_[i] = 0.f;
888  }
889 
890  // This works only for IP-related tracks
891 
892  flip_x_ = false;
893  flip_y_ = false;
894  float cotbeta = entry00_->cotbeta;
895  switch(Dtype_) {
896  case 0:
897  if(cotbeta < 0.f) {flip_y_ = true;}
898  break;
899  case 1:
900  if(locBz > 0.f) {
901  flip_y_ = true;
902  }
903  break;
904  case 2:
905  case 3:
906  case 4:
907  case 5:
908  if(locBx*locBz < 0.f) {
909  flip_x_ = true;
910  }
911  if(locBx < 0.f) {
912  flip_y_ = true;
913  }
914  break;
915  default:
916 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
917  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::illegal subdetector ID = " << iDtype << std::endl;
918 #else
919  std::cout << "SiPixelTemplate:2D:illegal subdetector ID = " << iDtype << std::endl;
920 #endif
921  }
922 
923  // Calculate signed quantities
924 
925  lorydrift_ = lorywidth_/2.;
926  if(flip_y_) lorydrift_ = -lorydrift_;
927  lorxdrift_ = lorxwidth_/2.;
928  if(flip_x_) lorxdrift_ = -lorxdrift_;
929 
930  for(int i=0; i<2 ; ++i) {
931  for(int j=0; j<5 ; ++j) {
932  // Charge loss switches sides when cot(beta) changes sign
933  if(flip_y_) {
934  xypary0x0_[1-i][j] = (float)entry00_->xypar[i][j];
935  xypary1x0_[1-i][j] = (float)entry00_->xypar[i][j];
936  xypary0x1_[1-i][j] = (float)entry00_->xypar[i][j];
937  lanpar_[1-i][j] = entry00_->lanpar[i][j];
938  } else {
939  xypary0x0_[i][j] = (float)entry00_->xypar[i][j];
940  xypary1x0_[i][j] = (float)entry00_->xypar[i][j];
941  xypary0x1_[i][j] = (float)entry00_->xypar[i][j];
942  lanpar_[i][j] = entry00_->lanpar[i][j];
943  }
944  }
945  }
946  return;
947 }
948 
949 
950 // *************************************************************************************************************************************
956 // *************************************************************************************************************************************
957 
958 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)
959 {
960  // Interpolate for a new set of track angles
961 
962  // Local variables
963  int pixx, pixy, k0, k1, l0, l1, deltax, deltay, iflipy, jflipx, imin, imax, jmin, jmax;
964  int m, n;
965  float dx, dy, ddx, ddy, adx, ady;
966  // const float deltaxy[2] = {8.33f, 12.5f};
967  const float deltaxy[2] = {16.67f, 25.0f};
968 
969  // Check to see if interpolation is valid
970 
971 
972  // next, determine the indices of the closest point in k (y-displacement), l (x-displacement)
973  // pixy and pixx are the indices of the struck pixel in the (Ty,Tx) system
974  // k0,k1 are the k-indices of the closest and next closest point
975  // l0,l1 are the l-indices of the closest and next closest point
976 
977  pixy = (int)floorf(yhit/ysize_);
978  dy = yhit-(pixy+0.5f)*ysize_;
979  if(flip_y_) {dy = -dy;}
980  k0 = (int)(dy/ysize_*6.f+3.5f);
981  if(k0 < 0) k0 = 0;
982  if(k0 > 6) k0 = 6;
983  ddy = 6.f*dy/ysize_ - (k0-3);
984  ady = fabs(ddy);
985  if(ddy > 0.f) {k1 = k0 + 1; if(k1 > 6) k1 = k0-1;} else {k1 = k0 - 1; if(k1 < 0) k1 = k0+1;}
986  pixx = (int)floorf(xhit/xsize_);
987  dx = xhit-(pixx+0.5f)*xsize_;
988  if(flip_x_) {dx = -dx;}
989  l0 = (int)(dx/xsize_*6.f+3.5f);
990  if(l0 < 0) l0 = 0;
991  if(l0 > 6) l0 = 6;
992  ddx = 6.f*dx/xsize_ - (l0-3);
993  adx = fabs(ddx);
994  if(ddx > 0.f) {l1 = l0 + 1; if(l1 > 6) l1 = l0-1;} else {l1 = l0 - 1; if(l1 < 0) l1 = l0+1;}
995 
996  // OK, lets do the template interpolation.
997 
998  // First find the limits of the indices for non-zero pixels
999 
1000  imin = std::min(entry00_->iymin,entry10_->iymin);
1001  imin_ = std::min(imin,entry01_->iymin);
1002 
1003  jmin = std::min(entry00_->jxmin,entry10_->jxmin);
1004  jmin_ = std::min(jmin,entry01_->jxmin);
1005 
1006  imax = std::max(entry00_->iymax,entry10_->iymax);
1007  imax_ = std::max(imax,entry01_->iymax);
1008 
1009  jmax = std::max(entry00_->jxmax,entry10_->jxmax);
1010  jmax_ = std::max(jmax,entry01_->jxmax);
1011 
1012  // Calculate the x and y offsets to make the new template
1013 
1014  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1015 
1016  ++pixy; ++pixx;
1017 
1018  // In the template store, the struck pixel is always (THy,THx)
1019 
1020  deltax = pixx - T2HX;
1021  deltay = pixy - T2HY;
1022 
1023  // First zero the local 2-d template
1024 
1025  for(int j=0; j<BXM2; ++j) {for(int i=0; i<BYM2; ++i) {xytemp_[j][i] = 0.f;}}
1026 
1027  // Loop over the non-zero part of the template index space and interpolate
1028 
1029  for(int j=jmin_; j<=jmax_; ++j) {
1030  // Flip indices as needed
1031  if(flip_x_) {jflipx=T2XSIZE-1-j; m = deltax+jflipx;} else {m = deltax+j;}
1032  for(int i=imin_; i<=imax_; ++i) {
1033  if(flip_y_) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
1034  if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
1035  xytemp_[m][n] = (float)entry00_->xytemp[k0][l0][i][j]
1036  + adx*(float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j])
1037  + ady*(float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1038  + adcota_*(float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1039  + adcotb_*(float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1040  }
1041  }
1042  }
1043 
1044  //combine rows and columns to simulate double pixels
1045 
1046  for(int n=1; n<BYM3; ++n) {
1047  if(ydouble[n]) {
1048  // Combine the y-columns
1049  for(int m=1; m<BXM3; ++m) {
1050  xytemp_[m][n] += xytemp_[m][n+1];
1051  }
1052  // Now shift the remaining pixels over by one column
1053  for(int i=n+1; i<BYM3; ++i) {
1054  for(int m=1; m<BXM3; ++m) {
1055  xytemp_[m][i] = xytemp_[m][i+1];
1056  }
1057  }
1058  }
1059  }
1060 
1061  //combine rows and columns to simulate double pixels
1062 
1063  for(int m=1; m<BXM3; ++m) {
1064  if(xdouble[m]) {
1065  // Combine the x-rows
1066  for(int n=1; n<BYM3; ++n) {
1067  xytemp_[m][n] += xytemp_[m+1][n];
1068  }
1069  // Now shift the remaining pixels over by one row
1070  for(int j=m+1; j<BXM3; ++j) {
1071  for(n=1; n<BYM3; ++n) {
1072  xytemp_[j][n] = xytemp_[j+1][n];
1073  }
1074  }
1075  }
1076  }
1077 
1078  // Finally, loop through and increment the external template
1079 
1080  float qtemptot = 0.f;
1081 
1082  for(int n=1; n<BYM3; ++n) {
1083  for(int m=1; m<BXM3; ++m) {
1084  if(xytemp_[m][n] != 0.f) {template2d[m][n] += xytemp_[m][n]; qtemptot += xytemp_[m][n];}
1085  }
1086  }
1087 
1088  QTemplate = qtemptot;
1089 
1090  if(derivatives) {
1091 
1092  float dxytempdx[2][BXM2][BYM2], dxytempdy[2][BXM2][BYM2];
1093 
1094  for(int k = 0; k<2; ++k) {
1095  for(int i = 0; i<BXM2; ++i) {
1096  for(int j = 0; j<BYM2; ++j) {
1097  dxytempdx[k][i][j] = 0.f;
1098  dxytempdy[k][i][j] = 0.f;
1099  dpdx2d[k][i][j] = 0.f;
1100  }
1101  }
1102  }
1103 
1104  // First do shifted +x template
1105 
1106  pixx = (int)floorf((xhit+deltaxy[0])/xsize_);
1107  dx = (xhit+deltaxy[0])-(pixx+0.5f)*xsize_;
1108  if(flip_x_) {dx = -dx;}
1109  l0 = (int)(dx/xsize_*6.f+3.5f);
1110  if(l0 < 0) l0 = 0;
1111  if(l0 > 6) l0 = 6;
1112  ddx = 6.f*dx/xsize_ - (l0-3);
1113  adx = fabs(ddx);
1114  if(ddx > 0.f) {l1 = l0 + 1; if(l1 > 6) l1 = l0-1;} else {l1 = l0 - 1; if(l1 < 0) l1 = l0+1;}
1115 
1116  // OK, lets do the template interpolation.
1117 
1118  // Calculate the x and y offsets to make the new template
1119 
1120  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1121 
1122  ++pixx;
1123 
1124  // In the template store, the struck pixel is always (THy,THx)
1125 
1126  deltax = pixx - T2HX;
1127 
1128  // Loop over the non-zero part of the template index space and interpolate
1129 
1130  for(int j=jmin_; j<=jmax_; ++j) {
1131  // Flip indices as needed
1132  if(flip_x_) {jflipx=T2XSIZE-1-j; m = deltax+jflipx;} else {m = deltax+j;}
1133  for(int i=imin_; i<=imax_; ++i) {
1134  if(flip_y_) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
1135  if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
1136  dxytempdx[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j]
1137  + adx*(float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j])
1138  + ady*(float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1139  + adcota_*(float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1140  + adcotb_*(float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1141  }
1142  }
1143  }
1144 
1145  //combine rows and columns to simulate double pixels
1146 
1147  for(int n=1; n<BYM3; ++n) {
1148  if(ydouble[n]) {
1149  // Combine the y-columns
1150  for(int m=1; m<BXM3; ++m) {
1151  dxytempdx[1][m][n] += dxytempdx[1][m][n+1];
1152  }
1153  // Now shift the remaining pixels over by one column
1154  for(int i=n+1; i<BYM3; ++i) {
1155  for(int m=1; m<BXM3; ++m) {
1156  dxytempdx[1][m][i] = dxytempdx[1][m][i+1];
1157  }
1158  }
1159  }
1160  }
1161 
1162  //combine rows and columns to simulate double pixels
1163 
1164  for(int m=1; m<BXM3; ++m) {
1165  if(xdouble[m]) {
1166  // Combine the x-rows
1167  for(int n=1; n<BYM3; ++n) {
1168  dxytempdx[1][m][n] += dxytempdx[1][m+1][n];
1169  }
1170  // Now shift the remaining pixels over by one row
1171  for(int j=m+1; j<BXM3; ++j) {
1172  for(int n=1; n<BYM3; ++n) {
1173  dxytempdx[1][j][n] = dxytempdx[1][j+1][n];
1174  }
1175  }
1176  }
1177  }
1178 
1179  // Next do shifted -x template
1180 
1181  pixx = (int)floorf((xhit-deltaxy[0])/xsize_);
1182  dx = (xhit-deltaxy[0])-(pixx+0.5f)*xsize_;
1183  if(flip_x_) {dx = -dx;}
1184  l0 = (int)(dx/xsize_*6.f+3.5f);
1185  if(l0 < 0) l0 = 0;
1186  if(l0 > 6) l0 = 6;
1187  ddx = 6.f*dx/xsize_ - (l0-3);
1188  adx = fabs(ddx);
1189  if(ddx > 0.f) {l1 = l0 + 1; if(l1 > 6) l1 = l0-1;} else {l1 = l0 - 1; if(l1 < 0) l1 = l0+1;}
1190 
1191  // OK, lets do the template interpolation.
1192 
1193  // Calculate the x and y offsets to make the new template
1194 
1195  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1196 
1197  ++pixx;
1198 
1199  // In the template store, the struck pixel is always (THy,THx)
1200 
1201  deltax = pixx - T2HX;
1202 
1203  // Loop over the non-zero part of the template index space and interpolate
1204 
1205  for(int j=jmin_; j<=jmax_; ++j) {
1206  // Flip indices as needed
1207  if(flip_x_) {jflipx=T2XSIZE-1-j; m = deltax+jflipx;} else {m = deltax+j;}
1208  for(int i=imin_; i<=imax_; ++i) {
1209  if(flip_y_) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
1210  if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
1211  dxytempdx[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j]
1212  + adx*(float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j])
1213  + ady*(float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1214  + adcota_*(float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1215  + adcotb_*(float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1216  }
1217  }
1218  }
1219 
1220  //combine rows and columns to simulate double pixels
1221 
1222  for(int n=1; n<BYM3; ++n) {
1223  if(ydouble[n]) {
1224  // Combine the y-columns
1225  for(int m=1; m<BXM3; ++m) {
1226  dxytempdx[0][m][n] += dxytempdx[0][m][n+1];
1227  }
1228  // Now shift the remaining pixels over by one column
1229  for(int i=n+1; i<BYM3; ++i) {
1230  for(int m=1; m<BXM3; ++m) {
1231  dxytempdx[0][m][i] = dxytempdx[0][m][i+1];
1232  }
1233  }
1234  }
1235  }
1236 
1237  //combine rows and columns to simulate double pixels
1238 
1239  for(int m=1; m<BXM3; ++m) {
1240  if(xdouble[m]) {
1241  // Combine the x-rows
1242  for(int n=1; n<BYM3; ++n) {
1243  dxytempdx[0][m][n] += dxytempdx[0][m+1][n];
1244  }
1245  // Now shift the remaining pixels over by one row
1246  for(int j=m+1; j<BXM3; ++j) {
1247  for(int n=1; n<BYM3; ++n) {
1248  dxytempdx[0][j][n] = dxytempdx[0][j+1][n];
1249  }
1250  }
1251  }
1252  }
1253 
1254 
1255  // Finally, normalize the derivatives and copy the results to the output array
1256 
1257  for(int n=1; n<BYM3; ++n) {
1258  for(int m=1; m<BXM3; ++m) {
1259  dpdx2d[0][m][n] = (dxytempdx[1][m][n] - dxytempdx[0][m][n])/(2.*deltaxy[0]);
1260  }
1261  }
1262 
1263  // Next, do shifted y template
1264 
1265  pixy = (int)floorf((yhit+deltaxy[1])/ysize_);
1266  dy = (yhit+deltaxy[1])-(pixy+0.5f)*ysize_;
1267  if(flip_y_) {dy = -dy;}
1268  k0 = (int)(dy/ysize_*6.f+3.5f);
1269  if(k0 < 0) k0 = 0;
1270  if(k0 > 6) k0 = 6;
1271  ddy = 6.f*dy/ysize_ - (k0-3);
1272  ady = fabs(ddy);
1273  if(ddy > 0.f) {k1 = k0 + 1; if(k1 > 6) k1 = k0-1;} else {k1 = k0 - 1; if(k1 < 0) k1 = k0+1;}
1274  pixx = (int)floorf(xhit/xsize_);
1275  dx = xhit-(pixx+0.5f)*xsize_;
1276  if(flip_x_) {dx = -dx;}
1277  l0 = (int)(dx/xsize_*6.f+3.5f);
1278  if(l0 < 0) l0 = 0;
1279  if(l0 > 6) l0 = 6;
1280  ddx = 6.f*dx/xsize_ - (l0-3);
1281  adx = fabs(ddx);
1282  if(ddx > 0.f) {l1 = l0 + 1; if(l1 > 6) l1 = l0-1;} else {l1 = l0 - 1; if(l1 < 0) l1 = l0+1;}
1283 
1284  // OK, lets do the template interpolation.
1285 
1286  // Calculate the x and y offsets to make the new template
1287 
1288  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1289 
1290  ++pixy; ++pixx;
1291 
1292  // In the template store, the struck pixel is always (THy,THx)
1293 
1294  deltax = pixx - T2HX;
1295  deltay = pixy - T2HY;
1296 
1297  // Loop over the non-zero part of the template index space and interpolate
1298 
1299  for(int j=jmin_; j<=jmax_; ++j) {
1300  // Flip indices as needed
1301  if(flip_x_) {jflipx=T2XSIZE-1-j; m = deltax+jflipx;} else {m = deltax+j;}
1302  for(int i=imin_; i<=imax_; ++i) {
1303  if(flip_y_) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
1304  if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
1305  dxytempdy[1][m][n] = (float)entry00_->xytemp[k0][l0][i][j]
1306  + adx*(float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j])
1307  + ady*(float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1308  + adcota_*(float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1309  + adcotb_*(float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1310  }
1311  }
1312  }
1313 
1314  //combine rows and columns to simulate double pixels
1315 
1316  for(int n=1; n<BYM3; ++n) {
1317  if(ydouble[n]) {
1318  // Combine the y-columns
1319  for(int m=1; m<BXM3; ++m) {
1320  dxytempdy[1][m][n] += dxytempdy[1][m][n+1];
1321  }
1322  // Now shift the remaining pixels over by one column
1323  for(int i=n+1; i<BYM3; ++i) {
1324  for(int m=1; m<BXM3; ++m) {
1325  dxytempdy[1][m][i] = dxytempdy[1][m][i+1];
1326  }
1327  }
1328  }
1329  }
1330 
1331  //combine rows and columns to simulate double pixels
1332 
1333  for(int m=1; m<BXM3; ++m) {
1334  if(xdouble[m]) {
1335  // Combine the x-rows
1336  for(int n=1; n<BYM3; ++n) {
1337  dxytempdy[1][m][n] += dxytempdy[1][m+1][n];
1338  }
1339  // Now shift the remaining pixels over by one row
1340  for(int j=m+1; j<BXM3; ++j) {
1341  for(int n=1; n<BYM3; ++n) {
1342  dxytempdy[1][j][n] = dxytempdy[1][j+1][n];
1343  }
1344  }
1345  }
1346  }
1347 
1348  // Next, do shifted -y template
1349 
1350  pixy = (int)floorf((yhit-deltaxy[1])/ysize_);
1351  dy = (yhit-deltaxy[1])-(pixy+0.5f)*ysize_;
1352  if(flip_y_) {dy = -dy;}
1353  k0 = (int)(dy/ysize_*6.f+3.5f);
1354  if(k0 < 0) k0 = 0;
1355  if(k0 > 6) k0 = 6;
1356  ddy = 6.f*dy/ysize_ - (k0-3);
1357  ady = fabs(ddy);
1358  if(ddy > 0.f) {k1 = k0 + 1; if(k1 > 6) k1 = k0-1;} else {k1 = k0 - 1; if(k1 < 0) k1 = k0+1;}
1359 
1360  // OK, lets do the template interpolation.
1361 
1362  // Calculate the x and y offsets to make the new template
1363 
1364  // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
1365 
1366  ++pixy;
1367 
1368  // In the template store, the struck pixel is always (THy,THx)
1369 
1370  deltay = pixy - T2HY;
1371 
1372  // Loop over the non-zero part of the template index space and interpolate
1373 
1374  for(int j=jmin_; j<=jmax_; ++j) {
1375  // Flip indices as needed
1376  if(flip_x_) {jflipx=T2XSIZE-1-j; m = deltax+jflipx;} else {m = deltax+j;}
1377  for(int i=imin_; i<=imax_; ++i) {
1378  if(flip_y_) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
1379  if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
1380  dxytempdy[0][m][n] = (float)entry00_->xytemp[k0][l0][i][j]
1381  + adx*(float)(entry00_->xytemp[k0][l1][i][j] - entry00_->xytemp[k0][l0][i][j])
1382  + ady*(float)(entry00_->xytemp[k1][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1383  + adcota_*(float)(entry01_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j])
1384  + adcotb_*(float)(entry10_->xytemp[k0][l0][i][j] - entry00_->xytemp[k0][l0][i][j]);
1385  }
1386  }
1387  }
1388 
1389  //combine rows and columns to simulate double pixels
1390 
1391  for(int n=1; n<BYM3; ++n) {
1392  if(ydouble[n]) {
1393  // Combine the y-columns
1394  for(int m=1; m<BXM3; ++m) {
1395  dxytempdy[0][m][n] += dxytempdy[0][m][n+1];
1396  }
1397  // Now shift the remaining pixels over by one column
1398  for(int i=n+1; i<BYM3; ++i) {
1399  for(int m=1; m<BXM3; ++m) {
1400  dxytempdy[0][m][i] = dxytempdy[0][m][i+1];
1401  }
1402  }
1403  }
1404  }
1405 
1406  //combine rows and columns to simulate double pixels
1407 
1408  for(int m=1; m<BXM3; ++m) {
1409  if(xdouble[m]) {
1410  // Combine the x-rows
1411  for(int n=1; n<BYM3; ++n) {
1412  dxytempdy[0][m][n] += dxytempdy[0][m+1][n];
1413  }
1414  // Now shift the remaining pixels over by one row
1415  for(int j=m+1; j<BXM3; ++j) {
1416  for(int n=1; n<BYM3; ++n) {
1417  dxytempdy[0][j][n] = dxytempdy[0][j+1][n];
1418  }
1419  }
1420  }
1421  }
1422 
1423  // Finally, normalize the derivatives and copy the results to the output array
1424 
1425  for(int n=1; n<BYM3; ++n) {
1426  for(int m=1; m<BXM3; ++m) {
1427  dpdx2d[1][m][n] = (dxytempdy[1][m][n] - dxytempdy[0][m][n])/(2.*deltaxy[1]);
1428  }
1429  }
1430  }
1431 
1432  return success_;
1433 } // xytemp
1434 
1435 
1436 
1437 // *************************************************************************************************************************************
1444 // *************************************************************************************************************************************
1445 
1446 bool SiPixelTemplate2D::xytemp(float xhit, float yhit, bool ydouble[BYM2], bool xdouble[BXM2], float template2d[BXM2][BYM2])
1447 {
1448  // Interpolate for a new set of track angles
1449 
1450  bool derivatives = false;
1451  float dpdx2d[2][BXM2][BYM2];
1452  float QTemplate;
1453 
1454  return SiPixelTemplate2D::xytemp(xhit, yhit, ydouble, xdouble, template2d, derivatives, dpdx2d, QTemplate);
1455 
1456 } // xytemp
1457 
1458 
1459 
1460 // *************************************************************************************************************************************
1470 // *************************************************************************************************************************************
1471 
1472 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])
1473 {
1474 
1475  // Local variables
1476 
1477  bool derivatives = false;
1478  float dpdx2d[2][BXM2][BYM2];
1479  float QTemplate;
1480  float locBx = 1.f;
1481  if(cotbeta < 0.f) {locBx = -1.f;}
1482  float locBz = locBx;
1483  if(cotalpha < 0.f) {locBz = -locBx;}
1484 
1485  bool yd[BYM2], xd[BXM2];
1486 
1487  yd[0] = false; yd[BYM2-1] = false;
1488  for(int i=0; i<TYSIZE; ++i) { yd[i+1] = ydouble[i];}
1489  xd[0] = false; xd[BXM2-1] = false;
1490  for(int j=0; j<TXSIZE; ++j) { xd[j+1] = xdouble[j];}
1491 
1492 
1493  // Interpolate for a new set of track angles
1494 
1495  if(SiPixelTemplate2D::interpolate(id, cotalpha, cotbeta, locBz, locBx)) {
1496  return SiPixelTemplate2D::xytemp(xhit, yhit, yd, xd, template2d, derivatives, dpdx2d, QTemplate);
1497  } else {
1498  return false;
1499  }
1500 
1501 } // xytemp
1502 
1503 
1504 // ************************************************************************************************************
1510 // ************************************************************************************************************
1511 void SiPixelTemplate2D::xysigma2(float qpixel, int index, float& xysig2)
1512 
1513 {
1514  // Interpolate using quantities already stored in the private variables
1515 
1516  // Local variables
1517  float sigi, sigi2, sigi3, sigi4, qscale, err2, err00;
1518 
1519  // Make sure that input is OK
1520 
1521 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1522  if(index < 1 || index >= BYM2) {
1523  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::ysigma2 called with index = " << index << std::endl;
1524  }
1525 #else
1526  assert(index > 0 && index < BYM2);
1527 #endif
1528 
1529  // Define the maximum signal to use in the parameterization
1530 
1531  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1532 
1533  if(qpixel < sxymax_) {
1534  sigi = qpixel;
1535  qscale = 1.f;
1536  } else {
1537  sigi = sxymax_;
1538  qscale = qpixel/sxymax_;
1539  }
1540  sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
1541  if(index <= T2HYP1) {
1542  err00 = xypary0x0_[0][0]+xypary0x0_[0][1]*sigi+xypary0x0_[0][2]*sigi2+xypary0x0_[0][3]*sigi3+xypary0x0_[0][4]*sigi4;
1543  err2 = err00
1544  +adcota_*(xypary0x1_[0][0]+xypary0x1_[0][1]*sigi+xypary0x1_[0][2]*sigi2+xypary0x1_[0][3]*sigi3+xypary0x1_[0][4]*sigi4 - err00)
1545  +adcotb_*(xypary1x0_[0][0]+xypary1x0_[0][1]*sigi+xypary1x0_[0][2]*sigi2+xypary1x0_[0][3]*sigi3+xypary1x0_[0][4]*sigi4 - err00);
1546  } else {
1547  err00 = xypary0x0_[1][0]+xypary0x0_[1][1]*sigi+xypary0x0_[1][2]*sigi2+xypary0x0_[1][3]*sigi3+xypary0x0_[1][4]*sigi4;
1548  err2 = err00
1549  +adcota_*(xypary0x1_[1][0]+xypary0x1_[1][1]*sigi+xypary0x1_[1][2]*sigi2+xypary0x1_[1][3]*sigi3+xypary0x1_[1][4]*sigi4 - err00)
1550  +adcotb_*(xypary1x0_[1][0]+xypary1x0_[1][1]*sigi+xypary1x0_[1][2]*sigi2+xypary1x0_[1][3]*sigi3+xypary1x0_[1][4]*sigi4 - err00);
1551  }
1552  xysig2 =qscale*err2;
1553  if(xysig2 <= 0.f) {xysig2 = s50_*s50_;}
1554 
1555  return;
1556 
1557 } // End xysigma2
1558 
1559 
1560 // ************************************************************************************************************
1562 // ************************************************************************************************************
1563 void SiPixelTemplate2D::landau_par(float lanpar[2][5])
1564 
1565 {
1566  // Interpolate using quantities already stored in the private variables
1567 
1568  for(int i=0; i<2; ++i) {
1569  for(int j=0; j<5; ++j) {
1570  lanpar[i][j] = lanpar_[i][j];
1571  }
1572  }
1573  return;
1574 
1575 } // End lan_par
1576 
1577 
1578 
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:277
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:958
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:1511
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
create_public_lumi_plots.in_file
in_file
Definition: create_public_lumi_plots.py:863
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:67
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:613
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:1563
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: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:841
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