CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiPixelGenError.cc
Go to the documentation of this file.
1 //
2 // SiPixelGenError.cc Version 1.00
3 //
4 // Object stores Lorentz widths, bias corrections, and errors for the Generic Algorithm
5 //
6 // Created by Morris Swartz on 10/27/06.
7 // Add some debugging messages. d.k. 5/14
8 //
9 
10 //#include <stdlib.h>
11 //#include <stdio.h>
12 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
13 #include <math.h>
14 #else
15 #include <math.h>
16 #endif
17 #include <algorithm>
18 #include <vector>
19 //#include "boost/multi_array.hpp"
20 #include <iostream>
21 #include <iomanip>
22 #include <sstream>
23 #include <fstream>
24 #include <list>
25 
26 
27 
28 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
29 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelGenError.h"
30 #include "FWCore/ParameterSet/interface/FileInPath.h"
31 #include "FWCore/MessageLogger/interface/MessageLogger.h"
32 #define LOGERROR(x) LogError(x)
33 #define LOGWARNING(x) LogWarning(x)
34 #define LOGINFO(x) LogInfo(x)
35 #define ENDL " "
36 #include "FWCore/Utilities/interface/Exception.h"
37 using namespace edm;
38 #else
39 #include "SiPixelGenError.h"
40 #define LOGERROR(x) std::cout << x << ": "
41 #define LOGINFO(x) std::cout << x << ": "
42 #define LOGWARNING(x) std::cout << x << ": "
43 #define ENDL std::endl
44 #endif
45 
46 //****************************************************************
51 //****************************************************************
52 bool SiPixelGenError::pushfile(int filenum, std::vector< SiPixelGenErrorStore > & thePixelTemp_)
53 {
54  // Add info stored in external file numbered filenum to theGenErrorStore
55 
56  // Local variables
57  int i, j, k;
58  float costrk[3]={0,0,0};
59  const char *tempfile;
60  // char title[80]; remove this
61  char c;
62  const int code_version={1};
63 
64 
65 
66  // Create a filename for this run
67 
68  std::ostringstream tout;
69 
70  // Create different path in CMSSW than standalone
71 
72 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
73  tout << "CalibTracker/SiPixelESProducers/data/generror_summary_zp"
74  << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
75  std::string tempf = tout.str();
76  edm::FileInPath file( tempf.c_str() );
77  tempfile = (file.fullPath()).c_str();
78 #else
79  tout << "generror_summary_zp" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
80  std::string tempf = tout.str();
81  tempfile = tempf.c_str();
82 #endif
83 
84  // open the Generic file
85 
86  std::ifstream in_file(tempfile, std::ios::in);
87 
88  if(in_file.is_open()) {
89 
90  // Create a local GenError storage entry
91 
92  SiPixelGenErrorStore theCurrentTemp;
93 
94  // Read-in a header string first and print it
95 
96  for (i=0; (c=in_file.get()) != '\n'; ++i) {
97  if(i < 79) {theCurrentTemp.head.title[i] = c;}
98  }
99  if(i > 78) {i=78;}
100  theCurrentTemp.head.title[i+1] ='\0';
101  LOGINFO("SiPixelGenError") << "Loading Pixel GenError File - " << theCurrentTemp.head.title << ENDL;
102 
103  // next, the header information
104 
105  in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >>
106  theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >>
107  theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
108  theCurrentTemp.head.zsize >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
109  theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
110 
111  if(in_file.fail()) {LOGERROR("SiPixelGenError") << "Error reading file, no GenError load" << ENDL; return false;}
112 
113  LOGINFO("SiPixelGenError") << "GenError ID = " << theCurrentTemp.head.ID << ", GenError Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
114  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
115  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
116  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
117  << ", 1/2 multi dcol threshold " << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold " << theCurrentTemp.head.ss50
118  << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", y Lorentz Bias " << theCurrentTemp.head.lorybias
119  << ", x Lorentz width " << theCurrentTemp.head.lorxwidth << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
120  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", " << theCurrentTemp.head.fbin[1]
121  << ", " << theCurrentTemp.head.fbin[2]
122  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
123 
124  if(theCurrentTemp.head.templ_version < code_version) {LOGERROR("SiPixelGenError") << "code expects version " << code_version << ", no GenError load" << ENDL; return false;}
125 
126 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
127 
128 // next, layout the 1-d/2-d structures needed to store GenError info
129 
130  theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
131 
132  theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
133 
134 #endif
135 
136 // next, loop over all y-angle entries
137 
138  for (i=0; i < theCurrentTemp.head.NTy; ++i) {
139 
140  in_file >> theCurrentTemp.enty[i].runnum >> costrk[0] >> costrk[1] >> costrk[2];
141 
142  if(in_file.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 1, no GenError load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
143 
144  // Calculate the alpha, beta, and cot(beta) for this entry
145 
146  theCurrentTemp.enty[i].cotalpha = costrk[0]/costrk[2];
147 
148  theCurrentTemp.enty[i].cotbeta = costrk[1]/costrk[2];
149 
150  in_file >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].dyone
151  >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
152 
153  if(in_file.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 2, no GenError load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
154 
155  in_file >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo
156  >> theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].qmin2;
157 
158  if(in_file.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 3, no GenError load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
159 
160  if(theCurrentTemp.enty[i].qmin <= 0.) {LOGERROR("SiPixelGenError") << "Error in GenError ID " << theCurrentTemp.head.ID << " qmin = " << theCurrentTemp.enty[i].qmin << ", run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
161 
162 
163  for (j=0; j<4; ++j) {
164 
165  in_file >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j];
166 
167  if(in_file.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 14a, no GenError load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
168  }
169 
170 
171  }
172 
173  // next, loop over all barrel x-angle entries
174 
175  for (k=0; k < theCurrentTemp.head.NTyx; ++k) {
176 
177  for (i=0; i < theCurrentTemp.head.NTxx; ++i) {
178 
179  in_file >> theCurrentTemp.entx[k][i].runnum >> costrk[0] >> costrk[1] >> costrk[2];
180 
181  if(in_file.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 17, no GenError load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
182 
183  // Calculate the alpha, beta, and cot(beta) for this entry
184 
185  theCurrentTemp.entx[k][i].cotalpha = costrk[0]/costrk[2];
186 
187  theCurrentTemp.entx[k][i].cotbeta = costrk[1]/costrk[2];
188 
189  in_file >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >> theCurrentTemp.entx[k][i].dyone
190  >> theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
191 
192  if(in_file.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 18, no GenError load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
193 
194  in_file >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >> theCurrentTemp.entx[k][i].dxtwo
195  >> theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].qmin2;
196  // >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav;
197 
198  if(in_file.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 19, no GenError load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
199 
200 
201  for (j=0; j<4; ++j) {
202 
203  in_file >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j];
204 
205  if(in_file.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 30a, no GenError load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
206  }
207 
208  }
209  }
210 
211 
212  in_file.close();
213 
214  // Add this info to the store
215 
216  thePixelTemp_.push_back(theCurrentTemp);
217 
218  postInit(thePixelTemp_);
219 
220  return true;
221 
222  } else {
223 
224  // If file didn't open, report this
225 
226  LOGERROR("SiPixelGenError") << "Error opening File" << tempfile << ENDL;
227  return false;
228 
229  }
230 
231 } // TempInit
232 
233 
234 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
235 
236 //****************************************************************
240 //****************************************************************
241 bool SiPixelGenError::pushfile(const SiPixelGenErrorDBObject& dbobject,
242  std::vector< SiPixelGenErrorStore > & thePixelTemp_) {
243  // Add GenError stored in external dbobject to theGenErrorStore
244 
245  // Local variables
246  int i, j, k;
247  float costrk[3]={0,0,0};
248  // const char *tempfile;
249  const int code_version={1};
250 
251  // We must create a new object because dbobject must be a const and our stream must not be
252  SiPixelGenErrorDBObject db = dbobject;
253 
254  // Create a local GenError storage entry
255  SiPixelGenErrorStore theCurrentTemp;
256 
257  // Fill the GenError storage for each GenError calibration stored in the db
258  for(int m=0; m<db.numOfTempl(); ++m) {
259 
260  // Read-in a header string first and print it
261 
262  SiPixelGenErrorDBObject::char2float temp;
263  for (i=0; i<20; ++i) {
264  temp.f = db.sVector()[db.index()];
265  theCurrentTemp.head.title[4*i] = temp.c[0];
266  theCurrentTemp.head.title[4*i+1] = temp.c[1];
267  theCurrentTemp.head.title[4*i+2] = temp.c[2];
268  theCurrentTemp.head.title[4*i+3] = temp.c[3];
269  db.incrementIndex(1);
270  }
271 
272  theCurrentTemp.head.title[79] = '\0';
273  LOGINFO("SiPixelGenError") << "Loading Pixel GenError File - "
274  << theCurrentTemp.head.title << ENDL;
275 
276  // next, the header information
277 
278  db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >>
279  theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale >> theCurrentTemp.head.s50 >>
280  theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >>
281  theCurrentTemp.head.zsize >> theCurrentTemp.head.ss50 >> theCurrentTemp.head.lorybias >> theCurrentTemp.head.lorxbias >>
282  theCurrentTemp.head.fbin[0] >> theCurrentTemp.head.fbin[1] >> theCurrentTemp.head.fbin[2];
283 
284  if(db.fail()) {LOGERROR("SiPixelGenError")
285  << "Error reading file, no GenError load" << ENDL; return false;}
286 
287  LOGINFO("SiPixelGenError") << "GenError ID = " << theCurrentTemp.head.ID << ", GenError Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
288  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
289  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
290  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
291  << ", 1/2 multi dcol threshold " << theCurrentTemp.head.s50 << ", 1/2 single dcol threshold " << theCurrentTemp.head.ss50
292  << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", y Lorentz Bias " << theCurrentTemp.head.lorybias
293  << ", x Lorentz width " << theCurrentTemp.head.lorxwidth << ", x Lorentz Bias " << theCurrentTemp.head.lorxbias
294  << ", Q/Q_avg fractions for Qbin defs " << theCurrentTemp.head.fbin[0] << ", " << theCurrentTemp.head.fbin[1]
295  << ", " << theCurrentTemp.head.fbin[2]
296  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
297 
298  LOGINFO("SiPixelGenError") << "Loading Pixel GenError - "
299  << theCurrentTemp.head.title << " version "
300  <<theCurrentTemp.head.templ_version <<" code v."
301  <<code_version<<ENDL;
302  if(theCurrentTemp.head.templ_version < code_version) {
303  LOGERROR("SiPixelGenError") << "code expects version " << code_version
304  << ", no GenError load" << ENDL; return false;}
305 
306 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
307 
308 // next, layout the 1-d/2-d structures needed to store GenError
309 
310  theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
311 
312  theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
313 
314 #endif
315 
316  // next, loop over all barrel y-angle entries
317 
318  for (i=0; i < theCurrentTemp.head.NTy; ++i) {
319 
320  db >> theCurrentTemp.enty[i].runnum >> costrk[0] >> costrk[1] >> costrk[2];
321 
322  if(db.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 1, no GenError load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
323 
324  // Calculate the alpha, beta, and cot(beta) for this entry
325 
326  theCurrentTemp.enty[i].cotalpha = costrk[0]/costrk[2];
327 
328  theCurrentTemp.enty[i].cotbeta = costrk[1]/costrk[2];
329 
330  db >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].pixmax >> theCurrentTemp.enty[i].dyone
331  >> theCurrentTemp.enty[i].syone >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone;
332 
333  if(db.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 2, no GenError load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
334 
335  db >> theCurrentTemp.enty[i].dytwo >> theCurrentTemp.enty[i].sytwo >> theCurrentTemp.enty[i].dxtwo
336  >> theCurrentTemp.enty[i].sxtwo >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].qmin2;
337 
338  for (j=0; j<4; ++j) {
339 
340  db >> theCurrentTemp.enty[i].yavggen[j] >> theCurrentTemp.enty[i].yrmsgen[j] >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j];
341 
342  if(db.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 14a, no GenError load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
343  }
344 
345  }
346 
347  // next, loop over all barrel x-angle entries
348 
349  for (k=0; k < theCurrentTemp.head.NTyx; ++k) {
350 
351  for (i=0; i < theCurrentTemp.head.NTxx; ++i) {
352 
353  db >> theCurrentTemp.entx[k][i].runnum >> costrk[0] >> costrk[1] >> costrk[2];
354 
355  if(db.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 17, no GenError load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
356 
357  // Calculate the alpha, beta, and cot(beta) for this entry
358 
359  theCurrentTemp.entx[k][i].cotalpha = costrk[0]/costrk[2];
360 
361  theCurrentTemp.entx[k][i].cotbeta = costrk[1]/costrk[2];
362 
363  db >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].pixmax >> theCurrentTemp.entx[k][i].dyone
364  >> theCurrentTemp.entx[k][i].syone >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone;
365 
366  if(db.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 18, no GenError load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
367 
368  db >> theCurrentTemp.entx[k][i].dytwo >> theCurrentTemp.entx[k][i].sytwo >> theCurrentTemp.entx[k][i].dxtwo
369  >> theCurrentTemp.entx[k][i].sxtwo >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].qmin2;
370 
371  for (j=0; j<4; ++j) {
372 
373  db >> theCurrentTemp.entx[k][i].yavggen[j] >> theCurrentTemp.entx[k][i].yrmsgen[j] >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j];
374 
375  if(db.fail()) {LOGERROR("SiPixelGenError") << "Error reading file 30a, no GenError load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
376  }
377 
378  }
379  }
380 
381 
382  // Add this GenError to the store
383 
384  thePixelTemp_.push_back(theCurrentTemp);
385 
386  postInit(thePixelTemp_);
387 
388  }
389  return true;
390 
391 } // TempInit
392 
393 #endif
394 
395 void SiPixelGenError::postInit(std::vector< SiPixelGenErrorStore > & thePixelTemp_) {
396 
397  for (auto & templ : thePixelTemp_) {
398  for ( auto iy=0; iy<templ.head.NTy; ++iy ) templ.cotbetaY[iy]=templ.enty[iy].cotbeta;
399  for ( auto iy=0; iy<templ.head.NTyx; ++iy ) templ.cotbetaX[iy]= templ.entx[iy][0].cotbeta;
400  for ( auto ix=0; ix<templ.head.NTxx; ++ix ) templ.cotalphaX[ix]=templ.entx[0][ix].cotalpha;
401  }
402 
403 }
404 
405 
406 
407 
408 // ************************************************************************************************************
430 // ************************************************************************************************************
431 // a simpler method just to return the LA
432 int SiPixelGenError::qbin(int id) {
433  // Find the index corresponding to id
434 
435  if(id != id_current_) {
436 
437  index_id_ = -1;
438  for( int i=0; i<(int)thePixelTemp_.size(); ++i) {
439  if(id == thePixelTemp_[i].head.ID) {
440  index_id_ = i;
441  id_current_ = id;
442  //
443  lorywidth_ = thePixelTemp_[i].head.lorywidth;
444  lorxwidth_ = thePixelTemp_[i].head.lorxwidth;
445  lorybias_ = thePixelTemp_[i].head.lorybias;
446  lorxbias_ = thePixelTemp_[i].head.lorxbias;
447 
448  //for(int j=0; j<3; ++j) {fbin_[j] = thePixelTemp_[i].head.fbin[j];}
449 
450  // Pixel sizes to the private variables
451  xsize_ = thePixelTemp_[i].head.xsize;
452  ysize_ = thePixelTemp_[i].head.ysize;
453  zsize_ = thePixelTemp_[i].head.zsize;
454 
455  break;
456  }
457  }
458  }
459  return index_id_;
460 }
461  //-----------------------------------------------------------------------
462  // Full method
463 int SiPixelGenError::qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus,
464  float& pixmx, float& sigmay, float& deltay, float& sigmax, float& deltax,
465  float& sy1, float& dy1, float& sy2, float& dy2, float& sx1, float& dx1,
466  float& sx2, float& dx2)
467 {
468  // Interpolate for a new set of track angles
469 
470 
471  // Find the index corresponding to id
472 
473 
474  if(id != id_current_) {
475 
476  index_id_ = -1;
477  for( int i=0; i<(int)thePixelTemp_.size(); ++i) {
478  if(id == thePixelTemp_[i].head.ID) {
479  index_id_ = i;
480  id_current_ = id;
481  lorywidth_ = thePixelTemp_[i].head.lorywidth;
482  lorxwidth_ = thePixelTemp_[i].head.lorxwidth;
483  lorybias_ = thePixelTemp_[i].head.lorybias;
484  lorxbias_ = thePixelTemp_[i].head.lorxbias;
485  for(int j=0; j<3; ++j) {fbin_[j] = thePixelTemp_[i].head.fbin[j];}
486 
487  // Pixel sizes to the private variables
488 
489  xsize_ = thePixelTemp_[i].head.xsize;
490  ysize_ = thePixelTemp_[i].head.ysize;
491  zsize_ = thePixelTemp_[i].head.zsize;
492 
493  break;
494  }
495  }
496  }
497 
498  int index = index_id_;
499 
500 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
501  if(index < 0 || index >= (int)thePixelTemp_.size()) {
502  throw cms::Exception("DataCorrupt") << "SiPixelGenError::qbin can't find needed GenError ID = " << id << std::endl;
503  }
504 #else
505  assert(index >= 0 && index < (int)thePixelTemp_.size());
506 #endif
507 
508  //
509 
510 
511  auto const & templ = thePixelTemp_[index];
512 
513  // Interpolate the absolute value of cot(beta)
514 
515  auto acotb = std::abs(cotbeta);
516 
517  // qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
518 
519  auto cotalpha0 = thePixelTemp_[index].enty[0].cotalpha;
520  auto qcorrect=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta*cotbeta+cotalpha0*cotalpha0));
521 
522  // for some cosmics, the ususal gymnastics are incorrect
523 
524  float cotb; bool flip_y;
525  if(thePixelTemp_[index].head.Dtype == 0) {
526  cotb = acotb;
527  flip_y = false;
528  if(cotbeta < 0.f) {flip_y = true;}
529  } else {
530  if(locBz < 0.f) {
531  cotb = cotbeta;
532  flip_y = false;
533  } else {
534  cotb = -cotbeta;
535  flip_y = true;
536  }
537  }
538 
539  // Copy the charge scaling factor to the private variable
540 
541  auto qscale = thePixelTemp_[index].head.qscale;
542 
543 
544  /*
545  lorywidth = thePixelTemp_[index].head.lorywidth;
546  if(locBz > 0.f) {lorywidth = -lorywidth;}
547  lorxwidth = thePixelTemp_[index].head.lorxwidth;
548  */
549 
550 
551  auto Ny = thePixelTemp_[index].head.NTy;
552  auto Nyx = thePixelTemp_[index].head.NTyx;
553  auto Nxx = thePixelTemp_[index].head.NTxx;
554 
555 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
556  if(Ny < 2 || Nyx < 1 || Nxx < 2) {
557  throw cms::Exception("DataCorrupt") << "GenError ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
558  }
559 #else
560  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
561 #endif
562 
563  // next, loop over all y-angle entries
564 
565  auto ilow = 0;
566  auto ihigh = 0;
567  auto yratio = 0.f;
568 
569  {
570  auto j = std::lower_bound(templ.cotbetaY,templ.cotbetaY+Ny,cotb);
571  if (j==templ.cotbetaY+Ny) { --j; yratio = 1.f; }
572  else if (j==templ.cotbetaY) { ++j; yratio = 0.f;}
573  else { yratio = (cotb - (*(j-1)) )/( (*j) - (*(j-1)) ) ; }
574 
575  ihigh = j-templ.cotbetaY;
576  ilow = ihigh-1;
577  }
578 
579 
580 
581  // Interpolate/store all y-related quantities (flip displacements when flip_y)
582 
583  dy1 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].dyone + yratio*thePixelTemp_[index].enty[ihigh].dyone;
584  if(flip_y) {dy1 = -dy1;}
585  sy1 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].syone + yratio*thePixelTemp_[index].enty[ihigh].syone;
586  dy2 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].dytwo + yratio*thePixelTemp_[index].enty[ihigh].dytwo;
587  if(flip_y) {dy2 = -dy2;}
588  sy2 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].sytwo + yratio*thePixelTemp_[index].enty[ihigh].sytwo;
589 
590  auto qavg = (1.f - yratio)*thePixelTemp_[index].enty[ilow].qavg + yratio*thePixelTemp_[index].enty[ihigh].qavg;
591  qavg *= qcorrect;
592  auto qmin = (1.f - yratio)*thePixelTemp_[index].enty[ilow].qmin + yratio*thePixelTemp_[index].enty[ihigh].qmin;
593  qmin *= qcorrect;
594  auto qmin2 = (1.f - yratio)*thePixelTemp_[index].enty[ilow].qmin2 + yratio*thePixelTemp_[index].enty[ihigh].qmin2;
595  qmin2 *= qcorrect;
596 
597 
598 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
599  if(qavg <= 0.f || qmin <= 0.f) {
600  throw cms::Exception("DataCorrupt") << "SiPixelGenError::qbin, qavg or qmin <= 0,"
601  << " Probably someone called the generic pixel reconstruction with an illegal trajectory state" << std::endl;
602  }
603 #else
604  assert(qavg > 0.f && qmin > 0.f);
605 #endif
606 
607  // Scale the input charge to account for differences between pixelav and CMSSW simulation or data
608 
609  auto qtotal = qscale*qclus;
610 
611  // uncertainty and final corrections depend upon total charge bin
612 
613  auto fq = qtotal/qavg;
614  int binq;
615  if(fq > fbin_[0]) {
616  binq=0;
617  } else {
618  if(fq > fbin_[1]) {
619  binq=1;
620  } else {
621  if(fq > fbin_[2]) {
622  binq=2;
623  } else {
624  binq=3;
625  }
626  }
627  }
628 
629  auto yavggen =(1.f - yratio)*thePixelTemp_[index].enty[ilow].yavggen[binq] + yratio*thePixelTemp_[index].enty[ihigh].yavggen[binq];
630  if(flip_y) {yavggen = -yavggen;}
631  auto yrmsgen =(1.f - yratio)*thePixelTemp_[index].enty[ilow].yrmsgen[binq] + yratio*thePixelTemp_[index].enty[ihigh].yrmsgen[binq];
632 
633 
634  // next, loop over all x-angle entries, first, find relevant y-slices
635 
636  auto iylow = 0;
637  auto iyhigh = 0;
638  auto yxratio = 0.f;
639 
640 
641  {
642  auto j = std::lower_bound(templ.cotbetaX,templ.cotbetaX+Nyx,acotb);
643  if (j==templ.cotbetaX+Nyx) { --j; yxratio = 1.f; }
644  else if (j==templ.cotbetaX) { ++j; yxratio = 0.f;}
645  else { yxratio = (acotb - (*(j-1)) )/( (*j) - (*(j-1)) ) ; }
646 
647  iyhigh = j-templ.cotbetaX;
648  iylow = iyhigh -1;
649  }
650 
651 
652 
653  ilow = ihigh = 0;
654  auto xxratio = 0.f;
655 
656  {
657  auto j = std::lower_bound(templ.cotalphaX,templ.cotalphaX+Nxx,cotalpha);
658  if (j==templ.cotalphaX+Nxx) { --j; xxratio = 1.f; }
659  else if (j==templ.cotalphaX) { ++j; xxratio = 0.f;}
660  else { xxratio = (cotalpha - (*(j-1)) )/( (*j) - (*(j-1)) ) ; }
661 
662  ihigh = j-templ.cotalphaX;
663  ilow = ihigh-1;
664  }
665 
666 
667 
668  dx1 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].dxone + xxratio*thePixelTemp_[index].entx[0][ihigh].dxone;
669  sx1 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].sxone + xxratio*thePixelTemp_[index].entx[0][ihigh].sxone;
670  dx2 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].dxtwo + xxratio*thePixelTemp_[index].entx[0][ihigh].dxtwo;
671  sx2 = (1.f - xxratio)*thePixelTemp_[index].entx[0][ilow].sxtwo + xxratio*thePixelTemp_[index].entx[0][ihigh].sxtwo;
672 
673  // pixmax is the maximum allowed pixel charge (used for truncation)
674 
675  pixmx=(1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].pixmax + xxratio*thePixelTemp_[index].entx[iylow][ihigh].pixmax)
676  +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].pixmax + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].pixmax);
677 
678  auto xavggen = (1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].xavggen[binq] + xxratio*thePixelTemp_[index].entx[iylow][ihigh].xavggen[binq])
679  +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].xavggen[binq] + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].xavggen[binq]);
680 
681  auto xrmsgen = (1.f - yxratio)*((1.f - xxratio)*thePixelTemp_[index].entx[iylow][ilow].xrmsgen[binq] + xxratio*thePixelTemp_[index].entx[iylow][ihigh].xrmsgen[binq])
682  +yxratio*((1.f - xxratio)*thePixelTemp_[index].entx[iyhigh][ilow].xrmsgen[binq] + xxratio*thePixelTemp_[index].entx[iyhigh][ihigh].xrmsgen[binq]);
683 
684 
685 
686  // Take the errors and bias from the correct charge bin
687 
688  sigmay = yrmsgen; deltay = yavggen;
689 
690  sigmax = xrmsgen; deltax = xavggen;
691 
692  // If the charge is too small (then flag it)
693 
694  if(qtotal < 0.95f*qmin) {binq = 5;} else {if(qtotal < 0.95f*qmin2) {binq = 4;}}
695 
696  return binq;
697 
698 } // qbin
i
int i
Definition: DBlmapReader.cc:9
SiPixelGenErrorEntry::qavg
float qavg
average cluster charge for this set of track angles (now includes threshold effects) ...
Definition: SiPixelGenError.h:28
SiPixelGenErrorHeader::ss50
float ss50
1/2 of the single hit dcol threshold in electrons
Definition: SiPixelGenError.h:65
SiPixelGenErrorEntry::qmin
float qmin
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
Definition: SiPixelGenError.h:38
SiPixelGenErrorHeader::Bfield
float Bfield
Bfield in Tesla.
Definition: SiPixelGenError.h:68
SiPixelGenErrorEntry::runnum
int runnum
&lt; Basic template entry corresponding to a single set of track angles
Definition: SiPixelGenError.h:25
SiPixelGenErrorEntry::yrmsgen
float yrmsgen[4]
generic algorithm: average y-rms of reconstruction binned in 4 charge bins
Definition: SiPixelGenError.h:41
SiPixelGenErrorEntry::dxone
float dxone
mean offset/correction for one pixel x-clusters
Definition: SiPixelGenError.h:32
LOGINFO
#define LOGINFO(x)
Definition: SiPixelGenError.cc:34
SiPixelGenErrorHeader::Vbias
float Vbias
detector bias potential in Volts
Definition: SiPixelGenError.h:61
SiPixelGenErrorEntry::xrmsgen
float xrmsgen[4]
generic algorithm: average x-rms of reconstruction binned in 4 charge bins
Definition: SiPixelGenError.h:43
assert
assert(m_qm.get())
EcalCondDB.db
tuple db
Definition: EcalCondDB.py:151
SiPixelGenErrorHeader::NTxx
int NTxx
number of Template x-entries in each slice
Definition: SiPixelGenError.h:54
SiPixelGenErrorHeader::fbin
float fbin[3]
The QBin definitions in Q_clus/Q_avg.
Definition: SiPixelGenError.h:69
SiPixelGenErrorEntry::syone
float syone
rms for one pixel y-clusters
Definition: SiPixelGenError.h:31
SiPixelGenErrorHeader::zsize
float zsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelGenError.h:72
SiPixelGenErrorStore::enty
SiPixelGenErrorEntry enty[60]
60 Barrel y templates spanning cluster lengths from 0px to +18px [28 entries for fpix] ...
Definition: SiPixelGenError.h:82
SiPixelGenErrorHeader::s50
float s50
1/2 of the multihit dcol threshold in electrons
Definition: SiPixelGenError.h:64
SiPixelGenErrorEntry::cotalpha
float cotalpha
cot(alpha) is proportional to cluster length in x and is basis of interpolation
Definition: SiPixelGenError.h:26
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
SiPixelGenErrorHeader::lorxbias
float lorxbias
estimate of x-lorentz bias
Definition: SiPixelGenError.h:60
create_public_lumi_plots.in_file
tuple in_file
Definition: create_public_lumi_plots.py:868
SiPixelGenErrorHeader::temperature
float temperature
detector temperature in deg K
Definition: SiPixelGenError.h:62
SiPixelGenErrorHeader::xsize
float xsize
pixel size (for future use in upgraded geometry)
Definition: SiPixelGenError.h:70
SiPixelGenError::postInit
static void postInit(std::vector< SiPixelGenErrorStore > &thePixelTemp_)
Definition: SiPixelGenError.cc:395
SiPixelGenErrorHeader::templ_version
int templ_version
Version number of the template to ensure code compatibility.
Definition: SiPixelGenError.h:67
SiPixelGenErrorEntry::cotbeta
float cotbeta
cot(beta) is proportional to cluster length in y and is basis of interpolation
Definition: SiPixelGenError.h:27
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
SiPixelGenErrorDBObject::sVector
std::vector< float > sVector() const
Definition: SiPixelGenErrorDBObject.h:67
SiPixelGenErrorHeader::ID
int ID
&lt; template header structure
Definition: SiPixelGenError.h:51
SiPixelGenErrorEntry::sxone
float sxone
rms for one pixel x-clusters
Definition: SiPixelGenError.h:33
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
j
int j
Definition: DBlmapReader.cc:9
SiPixelGenErrorHeader::fluence
float fluence
radiation fluence in n_eq/cm^2
Definition: SiPixelGenError.h:63
LOGERROR
#define LOGERROR(x)
Definition: SiPixelGenError.cc:32
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
SiPixelGenErrorEntry::dytwo
float dytwo
mean offset/correction for one double-pixel y-clusters
Definition: SiPixelGenError.h:34
SiPixelGenErrorHeader::NTy
int NTy
number of Template y entries
Definition: SiPixelGenError.h:52
SiPixelGenErrorHeader::qscale
float qscale
Charge scaling to match cmssw and pixelav.
Definition: SiPixelGenError.h:56
SiPixelGenErrorEntry::dxtwo
float dxtwo
mean offset/correction for one double-pixel x-clusters
Definition: SiPixelGenError.h:36
SiPixelGenErrorStore::head
SiPixelGenErrorHeader head
&lt; template storage structure
Definition: SiPixelGenError.h:78
SiPixelGenErrorHeader::ysize
float ysize
pixel size (for future use in upgraded geometry)
Definition: SiPixelGenError.h:71
SiPixelGenErrorHeader::title
char title[80]
template title
Definition: SiPixelGenError.h:66
SiPixelGenErrorEntry::sytwo
float sytwo
rms for one double-pixel y-clusters
Definition: SiPixelGenError.h:35
SiPixelGenErrorHeader::lorywidth
float lorywidth
estimate of y-lorentz width for optimal resolution
Definition: SiPixelGenError.h:57
ENDL
#define ENDL
Definition: SiPixelGenError.cc:35
SiPixelGenErrorEntry::yavggen
float yavggen[4]
generic algorithm: average y-bias of reconstruction binned in 4 charge bins
Definition: SiPixelGenError.h:40
SiPixelGenErrorEntry::dyone
float dyone
mean offset/correction for one pixel y-clusters
Definition: SiPixelGenError.h:30
SiPixelGenErrorHeader::Dtype
int Dtype
detector type (0=BPix, 1=FPix)
Definition: SiPixelGenError.h:55
SiPixelGenErrorHeader::NTyx
int NTyx
number of Template y-slices of x entries
Definition: SiPixelGenError.h:53
SiPixelGenErrorEntry::xavggen
float xavggen[4]
generic algorithm: average x-bias of reconstruction binned in 4 charge bins
Definition: SiPixelGenError.h:42
SiPixelGenError::pushfile
static bool pushfile(int filenum, std::vector< SiPixelGenErrorStore > &thePixelTemp_)
Definition: SiPixelGenError.cc:52
relval_steps.k
list k
Definition: relval_steps.py:1484
SiPixelGenErrorEntry::pixmax
float pixmax
maximum charge for individual pixels in cluster
Definition: SiPixelGenError.h:29
citk::for
for(const auto &isodef:isoDefs)
Definition: CITKPFIsolationSumProducer.cc:77
SiPixelGenErrorEntry::sxtwo
float sxtwo
rms for one double-pixel x-clusters
Definition: SiPixelGenError.h:37
SiPixelGenError::qbin
int qbin(int id, float cotalpha, float cotbeta, float locBz, float qclus, float &pixmx, float &sigmay, float &deltay, float &sigmax, float &deltax, float &sy1, float &dy1, float &sy2, float &dy2, float &sx1, float &dx1, float &sx2, float &dx2)
Definition: SiPixelGenError.cc:463
SiPixelGenErrorStore::entx
SiPixelGenErrorEntry entx[5][29]
29 Barrel x templates spanning cluster lengths from -6px (-1.125Rad) to +6px (+1.125Rad) in each of 5...
Definition: SiPixelGenError.h:83
SiPixelGenErrorHeader::lorybias
float lorybias
estimate of y-lorentz bias
Definition: SiPixelGenError.h:59
SiPixelGenErrorHeader::lorxwidth
float lorxwidth
estimate of x-lorentz width for optimal resolution
Definition: SiPixelGenError.h:58
Generated for CMSSW Reference Manual by doxygen 1.8.5