CMS 3D CMS Logo

SiPixelTemplate2D.cc
Go to the documentation of this file.
1 //
2 // SiPixelTemplate2D.cc Version 1.03
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 //
13 
14 //#include <stdlib.h>
15 //#include <stdio.h>
16 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
17 //#include <cmath.h>
18 #else
19 #include <math.h>
20 #endif
21 #include <algorithm>
22 #include <vector>
23 //#include "boost/multi_array.hpp"
24 #include <iostream>
25 #include <iomanip>
26 #include <sstream>
27 #include <fstream>
28 
29 
30 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
34 #define LOGERROR(x) LogError(x)
35 #define LOGINFO(x) LogInfo(x)
36 #define ENDL " "
38 using namespace edm;
39 #else
40 #include "SiPixelTemplate2D.h"
41 #define LOGERROR(x) std::cout << x << ": "
42 #define LOGINFO(x) std::cout << x << ": "
43 #define ENDL std::endl
44 #endif
45 
46 //****************************************************************
51 //****************************************************************
52 bool SiPixelTemplate2D::pushfile(int filenum, std::vector< SiPixelTemplateStore2D > & thePixelTemp_)
53 {
54  // Add template stored in external file numbered filenum to theTemplateStore
55 
56  // Local variables
57  int i, j, k, l, iy, jx;
58  const char *tempfile;
59  // char title[80]; remove this
60  char c;
61  const int code_version={16};
62 
63 
64 
65  // Create a filename for this run
66 
67  std::ostringstream tout;
68 
69  // Create different path in CMSSW than standalone
70 
71 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
72  tout << "CalibTracker/SiPixelESProducers/data/template_summary2D_zp"
73  << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
74  std::string tempf = tout.str();
75  edm::FileInPath file( tempf.c_str() );
76  tempfile = (file.fullPath()).c_str();
77 #else
78  tout << "template_summary2D_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
79  std::string tempf = tout.str();
80  tempfile = tempf.c_str();
81 #endif
82 
83  // open the template file
84 
85  std::ifstream in_file(tempfile, std::ios::in);
86 
87  if(in_file.is_open()) {
88 
89  // Create a local template storage entry
90 
91  SiPixelTemplateStore2D theCurrentTemp;
92 
93  // Read-in a header string first and print it
94 
95  for (i=0; (c=in_file.get()) != '\n'; ++i) {
96  if(i < 79) {theCurrentTemp.head.title[i] = c;}
97  }
98  if(i > 78) {i=78;}
99  theCurrentTemp.head.title[i+1] ='\0';
100  LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
101 
102  // next, the header information
103 
104  in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
105  >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
106  >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
107 
108  if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file, no template load" << ENDL; return false;}
109 
110  LOGINFO("SiPixelTemplate2D") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
111  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
112  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
113  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
114  << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
115  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
116 
117  if(theCurrentTemp.head.templ_version != code_version) {LOGERROR("SiPixelTemplate2D") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
118 
119  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;}
120 
121  // next, layout the 2-d structure needed to store template
122 
123  theCurrentTemp.entry.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
124 
125  // Read in the file info
126 
127  for (iy=0; iy < theCurrentTemp.head.NTyx; ++iy) {
128  for(jx=0; jx < theCurrentTemp.head.NTxx; ++jx) {
129 
130  in_file >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0]
131  >> theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2];
132 
133  if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
134 
135  // Calculate cot(alpha) and cot(beta) for this entry
136 
137  theCurrentTemp.entry[iy][jx].cotalpha = theCurrentTemp.entry[iy][jx].costrk[0]/theCurrentTemp.entry[iy][jx].costrk[2];
138 
139  theCurrentTemp.entry[iy][jx].cotbeta = theCurrentTemp.entry[iy][jx].costrk[1]/theCurrentTemp.entry[iy][jx].costrk[2];
140 
141  in_file >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >> theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin
142  >> theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >> theCurrentTemp.entry[iy][jx].jxmax;
143 
144  if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
145 
146  for (k=0; k<2; ++k) {
147 
148  in_file >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1]
149  >> theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >> theCurrentTemp.entry[iy][jx].xypar[k][4];
150 
151  if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
152 
153  }
154 
155  for (k=0; k<2; ++k) {
156 
157  in_file >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1]
158  >> theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >> theCurrentTemp.entry[iy][jx].lanpar[k][4];
159 
160  if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
161 
162  }
163 
164  for (l=0; l<7; ++l) {
165  for (k=0; k<7; ++k) {
166  for (j=0; j<T2XSIZE; ++j) {
167 
168  for (i=0; i<T2YSIZE; ++i) {in_file >> theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j];}
169 
170  if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
171  }
172  }
173  }
174 
175 
176  in_file >> theCurrentTemp.entry[iy][jx].chi2minone >> theCurrentTemp.entry[iy][jx].chi2avgone >> theCurrentTemp.entry[iy][jx].chi2min[0] >> theCurrentTemp.entry[iy][jx].chi2avg[0]
177  >> theCurrentTemp.entry[iy][jx].chi2min[1] >> theCurrentTemp.entry[iy][jx].chi2avg[1]>> theCurrentTemp.entry[iy][jx].chi2min[2] >> theCurrentTemp.entry[iy][jx].chi2avg[2]
178  >> theCurrentTemp.entry[iy][jx].chi2min[3] >> theCurrentTemp.entry[iy][jx].chi2avg[3];
179 
180  if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
181 
182  in_file >> theCurrentTemp.entry[iy][jx].spare[0] >> theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2] >> theCurrentTemp.entry[iy][jx].spare[3]
183  >> theCurrentTemp.entry[iy][jx].spare[4] >> theCurrentTemp.entry[iy][jx].spare[5] >> theCurrentTemp.entry[iy][jx].spare[6] >> theCurrentTemp.entry[iy][jx].spare[7]
184  >> theCurrentTemp.entry[iy][jx].spare[8] >> theCurrentTemp.entry[iy][jx].spare[9];
185 
186  if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
187 
188  in_file >> theCurrentTemp.entry[iy][jx].spare[10] >> theCurrentTemp.entry[iy][jx].spare[11] >> theCurrentTemp.entry[iy][jx].spare[12] >> theCurrentTemp.entry[iy][jx].spare[13]
189  >> theCurrentTemp.entry[iy][jx].spare[14] >> theCurrentTemp.entry[iy][jx].spare[15] >> theCurrentTemp.entry[iy][jx].spare[16] >> theCurrentTemp.entry[iy][jx].spare[17]
190  >> theCurrentTemp.entry[iy][jx].spare[18] >> theCurrentTemp.entry[iy][jx].spare[19];
191 
192  if(in_file.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
193  }
194 
195  }
196 
197 
198  in_file.close();
199 
200  // Add this template to the store
201 
202  thePixelTemp_.push_back(theCurrentTemp);
203 
204  return true;
205 
206  } else {
207 
208  // If file didn't open, report this
209 
210  LOGERROR("SiPixelTemplate2D") << "Error opening File" << tempfile << ENDL;
211  return false;
212 
213  }
214 
215 } // TempInit
216 
217 
218 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
219 
220 //****************************************************************
224 //****************************************************************
225 bool SiPixelTemplate2D::pushfile(const SiPixelTemplateDBObject& dbobject, std::vector< SiPixelTemplateStore2D > & thePixelTemp_)
226 {
227  // Add template stored in external dbobject to theTemplateStore
228 
229  // Local variables
230  int i, j, k, l, iy, jx;
231  // const char *tempfile;
232  const int code_version={16};
233 
234  // We must create a new object because dbobject must be a const and our stream must not be
235  SiPixelTemplateDBObject db = dbobject;
236 
237  // Create a local template storage entry
238  SiPixelTemplateStore2D theCurrentTemp;
239 
240  // Fill the template storage for each template calibration stored in the db
241  for(int m=0; m<db.numOfTempl(); ++m)
242  {
243 
244  // Read-in a header string first and print it
245 
247  for (i=0; i<20; ++i) {
248  temp.f = db.sVector()[db.index()];
249  theCurrentTemp.head.title[4*i] = temp.c[0];
250  theCurrentTemp.head.title[4*i+1] = temp.c[1];
251  theCurrentTemp.head.title[4*i+2] = temp.c[2];
252  theCurrentTemp.head.title[4*i+3] = temp.c[3];
253  db.incrementIndex(1);
254  }
255  theCurrentTemp.head.title[79] = '\0';
256  LOGINFO("SiPixelTemplate2D") << "Loading Pixel Template File - " << theCurrentTemp.head.title << ENDL;
257 
258  // next, the header information
259 
260  db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
261  >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
262  >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
263 
264  if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file, no template load" << ENDL; return false;}
265 
266  LOGINFO("SiPixelTemplate2D") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
267  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
268  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
269  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
270  << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
271  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
272 
273  if(theCurrentTemp.head.templ_version != code_version) {LOGERROR("SiPixelTemplate2D") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
274 
275  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;}
276 
277  // next, layout the 2-d structure needed to store template
278 
279  theCurrentTemp.entry.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
280 
281  // Read in the file info
282 
283  for (iy=0; iy < theCurrentTemp.head.NTyx; ++iy) {
284  for(jx=0; jx < theCurrentTemp.head.NTxx; ++jx) {
285 
286  db >> theCurrentTemp.entry[iy][jx].runnum >> theCurrentTemp.entry[iy][jx].costrk[0]
287  >> theCurrentTemp.entry[iy][jx].costrk[1] >> theCurrentTemp.entry[iy][jx].costrk[2];
288 
289  if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 1, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
290 
291  // Calculate cot(alpha) and cot(beta) for this entry
292 
293  theCurrentTemp.entry[iy][jx].cotalpha = theCurrentTemp.entry[iy][jx].costrk[0]/theCurrentTemp.entry[iy][jx].costrk[2];
294 
295  theCurrentTemp.entry[iy][jx].cotbeta = theCurrentTemp.entry[iy][jx].costrk[1]/theCurrentTemp.entry[iy][jx].costrk[2];
296 
297  db >> theCurrentTemp.entry[iy][jx].qavg >> theCurrentTemp.entry[iy][jx].pixmax >> theCurrentTemp.entry[iy][jx].sxymax >> theCurrentTemp.entry[iy][jx].iymin
298  >> theCurrentTemp.entry[iy][jx].iymax >> theCurrentTemp.entry[iy][jx].jxmin >> theCurrentTemp.entry[iy][jx].jxmax;
299 
300  if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 2, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
301 
302  for (k=0; k<2; ++k) {
303 
304  db >> theCurrentTemp.entry[iy][jx].xypar[k][0] >> theCurrentTemp.entry[iy][jx].xypar[k][1]
305  >> theCurrentTemp.entry[iy][jx].xypar[k][2] >> theCurrentTemp.entry[iy][jx].xypar[k][3] >> theCurrentTemp.entry[iy][jx].xypar[k][4];
306 
307  if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 3, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
308 
309  }
310 
311  for (k=0; k<2; ++k) {
312 
313  db >> theCurrentTemp.entry[iy][jx].lanpar[k][0] >> theCurrentTemp.entry[iy][jx].lanpar[k][1]
314  >> theCurrentTemp.entry[iy][jx].lanpar[k][2] >> theCurrentTemp.entry[iy][jx].lanpar[k][3] >> theCurrentTemp.entry[iy][jx].lanpar[k][4];
315 
316  if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 4, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
317 
318  }
319 
320  for (l=0; l<7; ++l) {
321  for (k=0; k<7; ++k) {
322  for (j=0; j<T2XSIZE; ++j) {
323 
324  for (i=0; i<T2YSIZE; ++i) {db >> theCurrentTemp.entry[iy][jx].xytemp[k][l][i][j];}
325 
326  if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 5, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
327  }
328  }
329  }
330 
331 
332  db >> theCurrentTemp.entry[iy][jx].chi2minone >> theCurrentTemp.entry[iy][jx].chi2avgone >> theCurrentTemp.entry[iy][jx].chi2min[0] >> theCurrentTemp.entry[iy][jx].chi2avg[0]
333  >> theCurrentTemp.entry[iy][jx].chi2min[1] >> theCurrentTemp.entry[iy][jx].chi2avg[1]>> theCurrentTemp.entry[iy][jx].chi2min[2] >> theCurrentTemp.entry[iy][jx].chi2avg[2]
334  >> theCurrentTemp.entry[iy][jx].chi2min[3] >> theCurrentTemp.entry[iy][jx].chi2avg[3];
335 
336  if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 6, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
337 
338  db >> theCurrentTemp.entry[iy][jx].spare[0] >> theCurrentTemp.entry[iy][jx].spare[1] >> theCurrentTemp.entry[iy][jx].spare[2] >> theCurrentTemp.entry[iy][jx].spare[3]
339  >> theCurrentTemp.entry[iy][jx].spare[4] >> theCurrentTemp.entry[iy][jx].spare[5] >> theCurrentTemp.entry[iy][jx].spare[6] >> theCurrentTemp.entry[iy][jx].spare[7]
340  >> theCurrentTemp.entry[iy][jx].spare[8] >> theCurrentTemp.entry[iy][jx].spare[9];
341 
342  if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 7, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
343 
344  db >> theCurrentTemp.entry[iy][jx].spare[10] >> theCurrentTemp.entry[iy][jx].spare[11] >> theCurrentTemp.entry[iy][jx].spare[12] >> theCurrentTemp.entry[iy][jx].spare[13]
345  >> theCurrentTemp.entry[iy][jx].spare[14] >> theCurrentTemp.entry[iy][jx].spare[15] >> theCurrentTemp.entry[iy][jx].spare[16] >> theCurrentTemp.entry[iy][jx].spare[17]
346  >> theCurrentTemp.entry[iy][jx].spare[18] >> theCurrentTemp.entry[iy][jx].spare[19];
347 
348  if(db.fail()) {LOGERROR("SiPixelTemplate2D") << "Error reading file 8, no template load, run # " << theCurrentTemp.entry[iy][jx].runnum << ENDL; return false;}
349 
350  }
351  }
352 
353  }
354 
355 
356 // Add this template to the store
357 
358  thePixelTemp_.push_back(theCurrentTemp);
359 
360  return true;
361 
362 } // TempInit
363 
364 #endif
365 
366 
367 
368 // *************************************************************************************************************************************
380 // *************************************************************************************************************************************
381 
382 bool SiPixelTemplate2D::xytemp(int id, float cotalpha, float cotbeta, float locBz, float xhit, float yhit, std::vector<bool>& ydouble, std::vector<bool>& xdouble, float template2d[BXM2][BYM2])
383 {
384  // Interpolate for a new set of track angles
385 
386  // Local variables
387  int i, j;
388  int pixx, pixy, k0, k1, l0, l1, imidx, deltax, deltay, iflipy, imin, imax, jmin, jmax;
389  int m, n;
390  float acotb, dcota, dcotb, dx, dy, ddx, ddy, adx, ady, tmpxy;
391  bool flip_y;
392 // std::vector <float> xrms(4), xgsig(4), xrmsc2m(4), xgsigc2m(4);
393  std::vector <float> chi2xavg(4), chi2xmin(4);
394 
395 
396 // Check to see if interpolation is valid
397 
398  if(id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
399 
400  cota_current_ = cotalpha; cotb_current_ = cotbeta; success_ = true;
401 
402  if(id != id_current_) {
403 
404  // Find the index corresponding to id
405 
406  index_id_ = -1;
407  for(i=0; i<(int)thePixelTemp_.size(); ++i) {
408 
409  if(id == thePixelTemp_[i].head.ID) {
410 
411  index_id_ = i;
412  id_current_ = id;
413 
414 // Copy the charge scaling factor to the private variable
415 
416  Dtype_ = thePixelTemp_[index_id_].head.Dtype;
417 
418 // Copy the charge scaling factor to the private variable
419 
420  qscale_ = thePixelTemp_[index_id_].head.qscale;
421 
422 // Copy the pseudopixel signal size to the private variable
423 
424  s50_ = thePixelTemp_[index_id_].head.s50;
425 
426 // Copy the Lorentz widths to private variables
427 
428  lorywidth_ = thePixelTemp_[index_id_].head.lorywidth;
429  lorxwidth_ = thePixelTemp_[index_id_].head.lorxwidth;
430 
431 // Copy the pixel sizes private variables
432 
433  xsize_ = thePixelTemp_[index_id_].head.xsize;
434  ysize_ = thePixelTemp_[index_id_].head.ysize;
435  zsize_ = thePixelTemp_[index_id_].head.zsize;
436 
437 // Determine the size of this template
438 
439  Nyx_ = thePixelTemp_[index_id_].head.NTyx;
440  Nxx_ = thePixelTemp_[index_id_].head.NTxx;
441 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
442  if(Nyx_ < 2 || Nxx_ < 2) {
443  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Nyx/Nxx = " << Nyx_ << "/" << Nxx_ << std::endl;
444  }
445 #else
446  assert(Nyx_ > 1 && Nxx_ > 1);
447 #endif
448  imidx = Nxx_/2;
449 
450  cotalpha0_ = thePixelTemp_[index_id_].entry[0][0].cotalpha;
451  cotalpha1_ = thePixelTemp_[index_id_].entry[0][Nxx_-1].cotalpha;
452  deltacota_ = (cotalpha1_-cotalpha0_)/(float)(Nxx_-1);
453 
454  cotbeta0_ = thePixelTemp_[index_id_].entry[0][imidx].cotbeta;
455  cotbeta1_ = thePixelTemp_[index_id_].entry[Nyx_-1][imidx].cotbeta;
456  deltacotb_ = (cotbeta1_-cotbeta0_)/(float)(Nyx_-1);
457 
458  break;
459  }
460  }
461  }
462  }
463 
464 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
465  if(index_id_ < 0 || index_id_ >= (int)thePixelTemp_.size()) {
466  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::interpolate can't find needed template ID = " << id
467  << ", Are you using the correct global tag?" << std::endl;
468  }
469 #else
470  assert(index_id_ >= 0 && index_id_ < (int)thePixelTemp_.size());
471 #endif
472 
473 // Check angle limits and et up interpolation parameters
474 
475  if(cotalpha < cotalpha0_) {
476  success_ = false;
477  jx0_ = 0;
478  jx1_ = 1;
479  adcota_ = 0.f;
480  } else if(cotalpha > cotalpha1_) {
481  success_ = false;
482  jx0_ = Nxx_ - 1;
483  jx1_ = jx0_ - 1;
484  adcota_ = 0.f;
485  } else {
486  jx0_ = (int)((cotalpha-cotalpha0_)/deltacota_+0.5f);
487  dcota = (cotalpha - (cotalpha0_ + jx0_*deltacota_))/deltacota_;
488  adcota_ = fabs(dcota);
489  if(dcota > 0.f) {jx1_ = jx0_ + 1;if(jx1_ > Nxx_-1) jx1_ = jx0_-1;} else {jx1_ = jx0_ - 1; if(jx1_ < 0) jx1_ = jx0_+1;}
490  }
491 
492 // Interpolate the absolute value of cot(beta)
493 
494  acotb = std::abs(cotbeta);
495 
496  if(acotb < cotbeta0_) {
497  success_ = false;
498  iy0_ = 0;
499  iy1_ = 1;
500  adcotb_ = 0.f;
501  } else if(acotb > cotbeta1_) {
502  success_ = false;
503  iy0_ = Nyx_ - 1;
504  iy1_ = iy0_ - 1;
505  adcotb_ = 0.f;
506  } else {
507  iy0_ = (int)((acotb-cotbeta0_)/deltacotb_+0.5f);
508  dcotb = (acotb - (cotbeta0_ + iy0_*deltacotb_))/deltacotb_;
509  adcotb_ = fabs(dcotb);
510  if(dcotb > 0.f) {iy1_ = iy0_ + 1; if(iy1_ > Nyx_-1) iy1_ = iy0_-1;} else {iy1_ = iy0_ - 1; if(iy1_ < 0) iy1_ = iy0_+1;}
511  }
512 
513 // This works only for IP-related tracks
514 
515  flip_y = false;
516  if(cotbeta < 0.f) {flip_y = true;}
517 
518 // If Fpix-related track has wrong cotbeta-field correlation, return false to trigger simple template instead
519 
520  if(Dtype_ == 1) {
521  if(cotbeta*locBz > 0.f) success_ = false;
522  }
523 
524 // Interpolate things in cot(alpha)-cot(beta)
525 
526  qavg_ = thePixelTemp_[index_id_].entry[iy0_][jx0_].qavg
527  +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].qavg - thePixelTemp_[index_id_].entry[iy0_][jx0_].qavg)
528  +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].qavg - thePixelTemp_[index_id_].entry[iy0_][jx0_].qavg);
529 
530  pixmax_ = thePixelTemp_[index_id_].entry[iy0_][jx0_].pixmax
531  +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].pixmax - thePixelTemp_[index_id_].entry[iy0_][jx0_].pixmax)
532  +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].pixmax - thePixelTemp_[index_id_].entry[iy0_][jx0_].pixmax);
533 
534  sxymax_ = thePixelTemp_[index_id_].entry[iy0_][jx0_].sxymax
535  +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].sxymax - thePixelTemp_[index_id_].entry[iy0_][jx0_].sxymax)
536  +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].sxymax - thePixelTemp_[index_id_].entry[iy0_][jx0_].sxymax);
537 
538  chi2avgone_ = thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avgone
539  +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].chi2avgone - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avgone)
540  +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].chi2avgone - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avgone);
541 
542  chi2minone_ = thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2minone
543  +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].chi2minone - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2minone)
544  +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].chi2minone - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2minone);
545 
546  for(i=0; i<4 ; ++i) {
547  chi2avg_[i] = thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avg[i]
548  +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].chi2avg[i] - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avg[i])
549  +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].chi2avg[i] - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2avg[i]);
550 
551  chi2min_[i] = thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2min[i]
552  +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].chi2min[i] - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2min[i])
553  +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].chi2min[i] - thePixelTemp_[index_id_].entry[iy0_][jx0_].chi2min[i]);
554  }
555 
556  for(i=0; i<2 ; ++i) {
557  for(j=0; j<5 ; ++j) {
558  // Charge loss switches sides when cot(beta) changes sign
559  if(flip_y) {
560  xypary0x0_[1-i][j] = thePixelTemp_[index_id_].entry[iy0_][jx0_].xypar[i][j];
561  xypary1x0_[1-i][j] = thePixelTemp_[index_id_].entry[iy1_][jx0_].xypar[i][j];
562  xypary0x1_[1-i][j] = thePixelTemp_[index_id_].entry[iy0_][jx1_].xypar[i][j];
563  lanpar_[1-i][j] = thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j]
564  +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].lanpar[i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j])
565  +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].lanpar[i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j]);
566  } else {
567  xypary0x0_[i][j] = thePixelTemp_[index_id_].entry[iy0_][jx0_].xypar[i][j];
568  xypary1x0_[i][j] = thePixelTemp_[index_id_].entry[iy1_][jx0_].xypar[i][j];
569  xypary0x1_[i][j] = thePixelTemp_[index_id_].entry[iy0_][jx1_].xypar[i][j];
570  lanpar_[i][j] = thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j]
571  +adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].lanpar[i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j])
572  +adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].lanpar[i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].lanpar[i][j]);
573  }
574  }
575  }
576 
577 // next, determine the indices of the closest point in k (y-displacement), l (x-displacement)
578 // pixy and pixx are the indices of the struck pixel in the (Ty,Tx) system
579 // k0,k1 are the k-indices of the closest and next closest point
580 // l0,l1 are the l-indices of the closest and next closest point
581 
582  pixy = (int)floorf(yhit/ysize_);
583  dy = yhit-(pixy+0.5f)*ysize_;
584  if(flip_y) {dy = -dy;}
585  k0 = (int)(dy/ysize_*6.f+3.5f);
586  if(k0 < 0) k0 = 0;
587  if(k0 > 6) k0 = 6;
588  ddy = 6.f*dy/ysize_ - (k0-3);
589  ady = fabs(ddy);
590  if(ddy > 0.f) {k1 = k0 + 1; if(k1 > 6) k1 = k0-1;} else {k1 = k0 - 1; if(k1 < 0) k1 = k0+1;}
591  pixx = (int)floorf(xhit/xsize_);
592  dx = xhit-(pixx+0.5f)*xsize_;
593  l0 = (int)(dx/xsize_*6.f+3.5f);
594  if(l0 < 0) l0 = 0;
595  if(l0 > 6) l0 = 6;
596  ddx = 6.f*dx/xsize_ - (l0-3);
597  adx = fabs(ddx);
598  if(ddx > 0.f) {l1 = l0 + 1; if(l1 > 6) l1 = l0-1;} else {l1 = l0 - 1; if(l1 < 0) l1 = l0+1;}
599 
600 // OK, lets do the template interpolation.
601 
602 // First find the limits of the indices for non-zero pixels
603 
604  imin = std::min(thePixelTemp_[index_id_].entry[iy0_][jx0_].iymin,thePixelTemp_[index_id_].entry[iy1_][jx0_].iymin);
605  imin = std::min(imin,thePixelTemp_[index_id_].entry[iy0_][jx1_].iymin);
606 
607  jmin = std::min(thePixelTemp_[index_id_].entry[iy0_][jx0_].jxmin,thePixelTemp_[index_id_].entry[iy1_][jx0_].jxmin);
608  jmin = std::min(jmin,thePixelTemp_[index_id_].entry[iy0_][jx1_].jxmin);
609 
610  imax = std::max(thePixelTemp_[index_id_].entry[iy0_][jx0_].iymax,thePixelTemp_[index_id_].entry[iy1_][jx0_].iymax);
611  imax = std::max(imax,thePixelTemp_[index_id_].entry[iy0_][jx1_].iymax);
612 
613  jmax = std::max(thePixelTemp_[index_id_].entry[iy0_][jx0_].jxmax,thePixelTemp_[index_id_].entry[iy1_][jx0_].jxmax);
614  jmax = std::max(jmax,thePixelTemp_[index_id_].entry[iy0_][jx1_].jxmax);
615 
616 // Calculate the x and y offsets to make the new template
617 
618 // First, shift the struck pixel coordinates to the (Ty+2, Tx+2) system
619 
620  ++pixy; ++pixx;
621 
622 // In the template store, the struck pixel is always (THy,THx)
623 
624  deltax = pixx - T2HX;
625  deltay = pixy - T2HY;
626 
627 // First zero the local 2-d template
628 
629  for(j=0; j<BXM2; ++j) {for(i=0; i<BYM2; ++i) {xytemp_[j][i] = 0.f;}}
630 
631 // Loop over the non-zero part of the template index space and interpolate
632 
633  for(j=jmin; j<=jmax; ++j) {
634  for(i=imin; i<=imax; ++i) {
635  m = deltax+j;
636 
637 // If cot(beta) < 0, we must flip the cluster, iflipy is the flipped y-index
638 
639  if(flip_y) {iflipy=T2YSIZE-1-i; n = deltay+iflipy;} else {n = deltay+i;}
640  if(m>=0 && m<=BXM3 && n>=0 && n<=BYM3) {
641  tmpxy = thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l0][i][j]
642  + adx*(thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l1][i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l0][i][j])
643  + ady*(thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k1][l0][i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l0][i][j])
644  + adcota_*(thePixelTemp_[index_id_].entry[iy0_][jx1_].xytemp[k0][l0][i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l0][i][j])
645  + adcotb_*(thePixelTemp_[index_id_].entry[iy1_][jx0_].xytemp[k0][l0][i][j] - thePixelTemp_[index_id_].entry[iy0_][jx0_].xytemp[k0][l0][i][j]);
646  if(tmpxy > 0.f) {xytemp_[m][n] = tmpxy;} else {xytemp_[m][n] = 0.f;}
647  }
648  }
649  }
650 
651 //combine rows and columns to simulate double pixels
652 
653  for(n=1; n<BYM3; ++n) {
654  if(ydouble[n-1]) {
655 // Combine the y-columns
656  for(m=1; m<BXM3; ++m) {
657  xytemp_[m][n] += xytemp_[m][n+1];
658  }
659 // Now shift the remaining pixels over by one column
660  for(i=n+1; i<BYM3; ++i) {
661  for(m=1; m<BXM3; ++m) {
662  xytemp_[m][i] = xytemp_[m][i+1];
663  }
664  }
665  }
666  }
667 
668 //combine rows and columns to simulate double pixels
669 
670  for(m=1; m<BXM3; ++m) {
671  if(xdouble[m-1]) {
672 // Combine the x-rows
673  for(n=1; n<BYM3; ++n) {
674  xytemp_[m][n] += xytemp_[m+1][n];
675  }
676 // Now shift the remaining pixels over by one row
677  for(j=m+1; j<BXM3; ++j) {
678  for(n=1; n<BYM3; ++n) {
679  xytemp_[j][n] = xytemp_[j+1][n];
680  }
681  }
682  }
683  }
684 
685 // Finally, loop through and increment the external template
686 
687  for(n=1; n<BYM3; ++n) {
688  for(m=1; m<BXM3; ++m) {
689  if(xytemp_[m][n] > 0.f) {template2d[m][n] += xytemp_[m][n];}
690  }
691  }
692 
693  return success_;
694 } // xytemp
695 
696 
697 
698 // *************************************************************************************************************************************
708 // *************************************************************************************************************************************
709 
710 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])
711 {
712  // Interpolate for a new set of track angles
713 
714  // Local variables
715  float locBz = -1;
716  if(cotbeta < 0.f) {locBz = -locBz;}
717 
718  return SiPixelTemplate2D::xytemp(id, cotalpha, cotbeta, locBz, xhit, yhit, ydouble, xdouble, template2d);
719 
720 } // xytemp
721 
722 
723 
724 // ************************************************************************************************************
730 // ************************************************************************************************************
731 void SiPixelTemplate2D::xysigma2(float qpixel, int index, float& xysig2)
732 
733 {
734  // Interpolate using quantities already stored in the private variables
735 
736  // Local variables
737  float sigi, sigi2, sigi3, sigi4, qscale, err2, err00;
738 
739  // Make sure that input is OK
740 
741 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
742  if(index < 2 || index >= BYM2) {
743  throw cms::Exception("DataCorrupt") << "SiPixelTemplate2D::ysigma2 called with index = " << index << std::endl;
744  }
745 #else
746  assert(index > 1 && index < BYM2);
747 #endif
748 
749  // Define the maximum signal to use in the parameterization
750 
751  // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
752 
753  if(qpixel < sxymax_) {
754  sigi = qpixel;
755  qscale = 1.f;
756  } else {
757  sigi = sxymax_;
758  qscale = qpixel/sxymax_;
759  }
760  sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
761  if(index <= THXP1) {
762  err00 = xypary0x0_[0][0]+xypary0x0_[0][1]*sigi+xypary0x0_[0][2]*sigi2+xypary0x0_[0][3]*sigi3+xypary0x0_[0][4]*sigi4;
763  err2 = err00
764  +adcota_*(xypary0x1_[0][0]+xypary0x1_[0][1]*sigi+xypary0x1_[0][2]*sigi2+xypary0x1_[0][3]*sigi3+xypary0x1_[0][4]*sigi4 - err00)
765  +adcotb_*(xypary1x0_[0][0]+xypary1x0_[0][1]*sigi+xypary1x0_[0][2]*sigi2+xypary1x0_[0][3]*sigi3+xypary1x0_[0][4]*sigi4 - err00);
766  } else {
767  err00 = xypary0x0_[1][0]+xypary0x0_[1][1]*sigi+xypary0x0_[1][2]*sigi2+xypary0x0_[1][3]*sigi3+xypary0x0_[1][4]*sigi4;
768  err2 = err00
769  +adcota_*(xypary0x1_[1][0]+xypary0x1_[1][1]*sigi+xypary0x1_[1][2]*sigi2+xypary0x1_[1][3]*sigi3+xypary0x1_[1][4]*sigi4 - err00)
770  +adcotb_*(xypary1x0_[1][0]+xypary1x0_[1][1]*sigi+xypary1x0_[1][2]*sigi2+xypary1x0_[1][3]*sigi3+xypary1x0_[1][4]*sigi4 - err00);
771  }
772  xysig2 =qscale*err2;
773  if(xysig2 <= 0.f) {LOGERROR("SiPixelTemplate2D") << "neg y-error-squared, id = " << id_current_ << ", index = " << index_id_ <<
774  ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_ << ", sigi = " << sigi << ENDL;}
775 
776  return;
777 
778 } // End xysigma2
779 
780 
781 // ************************************************************************************************************
783 // ************************************************************************************************************
784 void SiPixelTemplate2D::landau_par(float lanpar[2][5])
785 
786 {
787  // Interpolate using quantities already stored in the private variables
788 
789  // Local variables
790  int i,j;
791  for(i=0; i<2; ++i) {
792  for(j=0; j<5; ++j) {
793  lanpar[i][j] = lanpar_[i][j];
794  }
795  }
796  return;
797 
798 } // End lan_par
799 
800 
801 
float qscale
Charge scaling to match cmssw and pixelav.
int i
Definition: DBlmapReader.cc:9
float Vbias
detector bias potential in Volts
float lorxwidth
estimate of x-lorentz width from single pixel offset
float ysize
pixel size (for future use in upgraded geometry)
SiPixelTemplateHeader2D head
< template storage structure
void xysigma2(float qpixel, int index, float &xysig2)
#define T2HY
#define BXM3
#define T2HX
static bool pushfile(int filenum, std::vector< SiPixelTemplateStore2D > &thePixelTemp_)
#define ENDL
boost::multi_array< SiPixelTemplateEntry2D, 2 > entry
use 2d entry to store [47][5] barrel entries or [5][9] fpix entries
#define BXM2
#define T2XSIZE
int ID
template ID number
#define LOGERROR(x)
int NTxx
number of Template x-entries in each slice
#define LOGINFO(x)
float xsize
pixel size (for future use in upgraded geometry)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
#define BYM2
double f[11][100]
T min(T a, T b)
Definition: MathUtil.h:58
#define BYM3
char title[80]
< template header structure
int k[5][pyjets_maxn]
bool xytemp(int id, float cotalpha, float cotbeta, float locBz, float xhit, float yhit, std::vector< bool > &ydouble, std::vector< bool > &xdouble, float template2d[13+2][21+2])
float zsize
pixel size (for future use in upgraded geometry)
float lorywidth
estimate of y-lorentz width from single pixel offset
int NTy
number of Template y entries (= 0 for 2-D templates)
int Dtype
detector type (0=BPix, 1=FPix)
HLT enums.
void landau_par(float lanpar[2][5])
Return the Landau probability parameters for this set of cot(alpha, cot(beta)
std::vector< float > const & sVector() const
int templ_version
Version number of the template to ensure code compatibility.
float s50
1/2 of the readout threshold in ADC units
float temperature
detector temperature in deg K
float Bfield
Bfield in Tesla.
#define T2YSIZE
float fluence
radiation fluence in n_eq/cm^2
int NTyx
number of Template y-slices of x entries
#define THXP1