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