CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiStripTemplate.cc
Go to the documentation of this file.
1 //
2 // SiStripTemplate.cc Version 1.00 (based on SiPixelTemplate v8.20)
3 //
4 // V1.05 - add VI optimizations from pixel template object
5 // V1.06 - increase angular acceptance (and structure size)
6 // V2.00 - add barycenter interpolation and getters, fix calculation for charge deposition to accommodate cota-offsets in the central cotb entries.
7 // V2.01 - fix problem with number of spare entries
8 //
9 
10 // Created by Morris Swartz on 2/2/11.
11 //
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 #include <list>
29 
30 
31 
32 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
33 #include "RecoLocalTracker/SiStripRecHitConverter/interface/SiStripTemplate.h"
34 #include "FWCore/ParameterSet/interface/FileInPath.h"
35 #include "FWCore/MessageLogger/interface/MessageLogger.h"
36 #define LOGERROR(x) LogError(x)
37 #define LOGINFO(x) LogInfo(x)
38 #define ENDL " "
39 #include "FWCore/Utilities/interface/Exception.h"
40 using namespace edm;
41 #else
42 #include "SiStripTemplate.h"
43 #define LOGERROR(x) std::cout << x << ": "
44 #define LOGINFO(x) std::cout << x << ": "
45 #define ENDL std::endl
46 #endif
47 
48 //****************************************************************
53 //****************************************************************
54 bool SiStripTemplate::pushfile(int filenum)
55 {
56  // Add template stored in external file numbered filenum to theTemplateStore
57 
58  // Local variables
59  int i, j, k, l;
60  float qavg_avg;
61  const char *tempfile;
62  // char title[80]; remove this
63  char c;
64  const int code_version={18};
65 
66 
67 
68  // Create a filename for this run
69 
70  std::ostringstream tout;
71 
72  // Create different path in CMSSW than standalone
73 
74 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
75  tout << "CalibTracker/SiPixelESProducers/data/stemplate_summary_p" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
76  std::string tempf = tout.str();
77  edm::FileInPath file( tempf.c_str() );
78  tempfile = (file.fullPath()).c_str();
79 #else
80  tout << "stemplate_summary_p" << std::setw(4) << std::setfill('0') << std::right << filenum << ".out" << std::ends;
81  std::string tempf = tout.str();
82  tempfile = tempf.c_str();
83 #endif
84 
85  // open the template file
86 
87  std::ifstream in_file(tempfile, std::ios::in);
88 
89  if(in_file.is_open()) {
90 
91  // Create a local template storage entry
92 
93  SiStripTemplateStore theCurrentTemp;
94 
95  // Read-in a header string first and print it
96 
97  for (i=0; (c=in_file.get()) != '\n'; ++i) {
98  if(i < 79) {theCurrentTemp.head.title[i] = c;}
99  }
100  if(i > 78) {i=78;}
101  theCurrentTemp.head.title[i+1] ='\0';
102  LOGINFO("SiStripTemplate") << "Loading Strip Template File - " << theCurrentTemp.head.title << ENDL;
103 
104  // next, the header information
105 
106  in_file >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
107  >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
108  >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
109 
110  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file, no template load" << ENDL; return false;}
111 
112  LOGINFO("SiStripTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
113  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
114  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
115  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
116  << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
117  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
118 
119  if(theCurrentTemp.head.templ_version < code_version) {LOGERROR("SiStripTemplate") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
120 
121 #ifdef SI_STRIP_TEMPLATE_USE_BOOST
122 
123 // next, layout the 1-d/2-d structures needed to store template
124 
125  theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
126 
127  theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
128 
129 #endif
130 
131 // next, loop over all y-angle entries
132 
133  for (i=0; i < theCurrentTemp.head.NTy; ++i) {
134 
135  in_file >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0]
136  >> theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2];
137 
138  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
139 
140  // Calculate the alpha, beta, and cot(beta) for this entry
141 
142  theCurrentTemp.enty[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
143 
144  theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0]/theCurrentTemp.enty[i].costrk[2];
145 
146  theCurrentTemp.enty[i].beta = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
147 
148  theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1]/theCurrentTemp.enty[i].costrk[2];
149 
150  in_file >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].sxmax >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clslenx;
151 
152  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
153  for (j=0; j<2; ++j) {
154 
155  in_file >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1]
156  >> theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
157 
158  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 6, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
159 
160  }
161 
162  qavg_avg = 0.f;
163  for (j=0; j<9; ++j) {
164 
165  for (k=0; k<TSXSIZE; ++k) {in_file >> theCurrentTemp.enty[i].xtemp[j][k]; qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];}
166 
167  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 7, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
168  }
169  theCurrentTemp.enty[i].qavg_avg = qavg_avg/9.;
170 
171  for (j=0; j<4; ++j) {
172 
173  in_file >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >> theCurrentTemp.enty[i].xgsig[j];
174 
175  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 10, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
176  }
177 
178  for (j=0; j<4; ++j) {
179 
180  in_file >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >> theCurrentTemp.enty[i].xflpar[j][2]
181  >> theCurrentTemp.enty[i].xflpar[j][3] >> theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
182 
183  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 11, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
184  }
185 
186  for (j=0; j<4; ++j) {
187 
188  in_file >> theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j] >> theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
189 
190  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 12, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
191  }
192 
193  for (j=0; j<4; ++j) {
194 
195  in_file >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >> theCurrentTemp.enty[i].xgx0c2m[j] >> theCurrentTemp.enty[i].xgsigc2m[j];
196 
197  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
198  }
199 
200  for (j=0; j<4; ++j) {
201 
202  in_file >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >> theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
203 
204  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14b, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
205  }
206 
207  for (j=0; j<4; ++j) {
208 
209  in_file >> theCurrentTemp.enty[i].xavgbcn[j] >> theCurrentTemp.enty[i].xrmsbcn[j] >> theCurrentTemp.enty[i].xgx0bcn[j] >> theCurrentTemp.enty[i].xgsigbcn[j];
210 
211  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14c, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
212  }
213 
214  in_file >> theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2
215  >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav
216  >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >> theCurrentTemp.enty[i].kappavav2 >> theCurrentTemp.enty[i].spare[0];
217 
218  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 15, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
219 
220  in_file >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1] >> theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracxone
221  >> theCurrentTemp.enty[i].spare[1] >> theCurrentTemp.enty[i].spare[2] >> theCurrentTemp.enty[i].spare[3] >> theCurrentTemp.enty[i].spare[4]
222  >> theCurrentTemp.enty[i].spare[5] >> theCurrentTemp.enty[i].spare[6];
223 
224  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 16, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
225 
226  }
227 
228  // next, loop over all barrel x-angle entries
229 
230  for (k=0; k < theCurrentTemp.head.NTyx; ++k) {
231 
232  for (i=0; i < theCurrentTemp.head.NTxx; ++i) {
233 
234  in_file >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0]
235  >> theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2];
236 
237  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 17, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
238 
239  // Calculate the alpha, beta, and cot(beta) for this entry
240 
241  theCurrentTemp.entx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
242 
243  theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0]/theCurrentTemp.entx[k][i].costrk[2];
244 
245  theCurrentTemp.entx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
246 
247  theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1]/theCurrentTemp.entx[k][i].costrk[2];
248 
249  in_file >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clslenx;
250 
251  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 18, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
252 
253  for (j=0; j<2; ++j) {
254 
255  in_file >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1]
256  >> theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >> theCurrentTemp.entx[k][i].xpar[j][4];
257 
258  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 19, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
259 
260  }
261 
262  qavg_avg = 0.f;
263  for (j=0; j<9; ++j) {
264 
265  for (l=0; l<TSXSIZE; ++l) {in_file >> theCurrentTemp.entx[k][i].xtemp[j][l]; qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];}
266 
267  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 20, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
268  }
269  theCurrentTemp.entx[k][i].qavg_avg = qavg_avg/9.;
270 
271  for (j=0; j<4; ++j) {
272 
273  in_file >> theCurrentTemp.entx[k][i].xavg[j] >> theCurrentTemp.entx[k][i].xrms[j] >> theCurrentTemp.entx[k][i].xgx0[j] >> theCurrentTemp.entx[k][i].xgsig[j];
274 
275  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 21, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
276  }
277 
278  for (j=0; j<4; ++j) {
279 
280  in_file >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >> theCurrentTemp.entx[k][i].xflpar[j][2]
281  >> theCurrentTemp.entx[k][i].xflpar[j][3] >> theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
282 
283  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 22, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
284  }
285 
286  for (j=0; j<4; ++j) {
287 
288  in_file >> theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j] >> theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
289 
290  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 23, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
291  }
292 
293  for (j=0; j<4; ++j) {
294 
295  in_file >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >> theCurrentTemp.entx[k][i].xgx0c2m[j] >> theCurrentTemp.entx[k][i].xgsigc2m[j];
296 
297  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 24, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
298  }
299 
300  for (j=0; j<4; ++j) {
301 
302  in_file >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j] >> theCurrentTemp.entx[k][i].xgx0gen[j] >> theCurrentTemp.entx[k][i].xgsiggen[j];
303 
304  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 25, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
305  }
306 
307  for (j=0; j<4; ++j) {
308 
309  in_file >> theCurrentTemp.entx[k][i].xavgbcn[j] >> theCurrentTemp.entx[k][i].xrmsbcn[j] >> theCurrentTemp.entx[k][i].xgx0bcn[j] >> theCurrentTemp.entx[k][i].xgsigbcn[j];
310 
311  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 26, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
312  }
313 
314  in_file >> theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >> theCurrentTemp.entx[k][i].qmin2
315  >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav
316  >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >> theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].spare[0];
317 
318  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 27, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
319 
320  in_file >> theCurrentTemp.entx[k][i].qbfrac[0] >> theCurrentTemp.entx[k][i].qbfrac[1] >> theCurrentTemp.entx[k][i].qbfrac[2] >> theCurrentTemp.entx[k][i].fracxone
321  >> theCurrentTemp.entx[k][i].spare[1] >> theCurrentTemp.entx[k][i].spare[2] >> theCurrentTemp.entx[k][i].spare[3] >> theCurrentTemp.entx[k][i].spare[4]
322  >> theCurrentTemp.entx[k][i].spare[5] >> theCurrentTemp.entx[k][i].spare[6];
323 
324  if(in_file.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 28, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
325 
326 
327  }
328  }
329 
330 
331  in_file.close();
332 
333  // Add this template to the store
334 
335  theStripTemp_.push_back(theCurrentTemp);
336 
337  return true;
338 
339  } else {
340 
341  // If file didn't open, report this
342 
343  LOGERROR("SiStripTemplate") << "Error opening File " << tempfile << ENDL;
344  return false;
345 
346  }
347 
348 } // TempInit
349 
350 
351 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
352 
353 //****************************************************************
357 //****************************************************************
358 bool SiStripTemplate::pushfile(const SiPixelTemplateDBObject& dbobject)
359 {
360  // Add template stored in external dbobject to theTemplateStore
361 
362  // Local variables
363  int i, j, k, l;
364  float qavg_avg;
365  // const char *tempfile;
366  const int code_version={17};
367 
368  // We must create a new object because dbobject must be a const and our stream must not be
369  SiPixelTemplateDBObject db = dbobject;
370 
371  // Create a local template storage entry
372  SiStripTemplateStore theCurrentTemp;
373 
374  // Fill the template storage for each template calibration stored in the db
375  for(int m=0; m<db.numOfTempl(); ++m)
376  {
377 
378  // Read-in a header string first and print it
379 
380  SiPixelTemplateDBObject::char2float temp;
381  for (i=0; i<20; ++i) {
382  temp.f = db.sVector()[db.index()];
383  theCurrentTemp.head.title[4*i] = temp.c[0];
384  theCurrentTemp.head.title[4*i+1] = temp.c[1];
385  theCurrentTemp.head.title[4*i+2] = temp.c[2];
386  theCurrentTemp.head.title[4*i+3] = temp.c[3];
387  db.incrementIndex(1);
388  }
389  theCurrentTemp.head.title[79] = '\0';
390  LOGINFO("SiStripTemplate") << "Loading Strip Template File - " << theCurrentTemp.head.title << ENDL;
391 
392  // next, the header information
393 
394  db >> theCurrentTemp.head.ID >> theCurrentTemp.head.templ_version >> theCurrentTemp.head.Bfield >> theCurrentTemp.head.NTy >> theCurrentTemp.head.NTyx >> theCurrentTemp.head.NTxx
395  >> theCurrentTemp.head.Dtype >> theCurrentTemp.head.Vbias >> theCurrentTemp.head.temperature >> theCurrentTemp.head.fluence >> theCurrentTemp.head.qscale
396  >> theCurrentTemp.head.s50 >> theCurrentTemp.head.lorywidth >> theCurrentTemp.head.lorxwidth >> theCurrentTemp.head.ysize >> theCurrentTemp.head.xsize >> theCurrentTemp.head.zsize;
397 
398  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file, no template load" << ENDL; return false;}
399 
400  LOGINFO("SiStripTemplate") << "Template ID = " << theCurrentTemp.head.ID << ", Template Version " << theCurrentTemp.head.templ_version << ", Bfield = " << theCurrentTemp.head.Bfield
401  << ", NTy = " << theCurrentTemp.head.NTy << ", NTyx = " << theCurrentTemp.head.NTyx<< ", NTxx = " << theCurrentTemp.head.NTxx << ", Dtype = " << theCurrentTemp.head.Dtype
402  << ", Bias voltage " << theCurrentTemp.head.Vbias << ", temperature "
403  << theCurrentTemp.head.temperature << ", fluence " << theCurrentTemp.head.fluence << ", Q-scaling factor " << theCurrentTemp.head.qscale
404  << ", 1/2 threshold " << theCurrentTemp.head.s50 << ", y Lorentz Width " << theCurrentTemp.head.lorywidth << ", x Lorentz width " << theCurrentTemp.head.lorxwidth
405  << ", pixel x-size " << theCurrentTemp.head.xsize << ", y-size " << theCurrentTemp.head.ysize << ", zsize " << theCurrentTemp.head.zsize << ENDL;
406 
407  if(theCurrentTemp.head.templ_version < code_version) {LOGERROR("SiStripTemplate") << "code expects version " << code_version << ", no template load" << ENDL; return false;}
408 
409 
410 #ifdef SI_PIXEL_TEMPLATE_USE_BOOST
411 
412 // next, layout the 1-d/2-d structures needed to store template
413 
414  theCurrentTemp.enty.resize(boost::extents[theCurrentTemp.head.NTy]);
415 
416  theCurrentTemp.entx.resize(boost::extents[theCurrentTemp.head.NTyx][theCurrentTemp.head.NTxx]);
417 
418 #endif
419 
420  // next, loop over all y-angle entries
421 
422  for (i=0; i < theCurrentTemp.head.NTy; ++i) {
423 
424  db >> theCurrentTemp.enty[i].runnum >> theCurrentTemp.enty[i].costrk[0]
425  >> theCurrentTemp.enty[i].costrk[1] >> theCurrentTemp.enty[i].costrk[2];
426 
427  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 1, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
428 
429  // Calculate the alpha, beta, and cot(beta) for this entry
430 
431  theCurrentTemp.enty[i].alpha = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[0]));
432 
433  theCurrentTemp.enty[i].cotalpha = theCurrentTemp.enty[i].costrk[0]/theCurrentTemp.enty[i].costrk[2];
434 
435  theCurrentTemp.enty[i].beta = static_cast<float>(atan2((double)theCurrentTemp.enty[i].costrk[2], (double)theCurrentTemp.enty[i].costrk[1]));
436 
437  theCurrentTemp.enty[i].cotbeta = theCurrentTemp.enty[i].costrk[1]/theCurrentTemp.enty[i].costrk[2];
438 
439  db >> theCurrentTemp.enty[i].qavg >> theCurrentTemp.enty[i].sxmax >> theCurrentTemp.enty[i].dxone >> theCurrentTemp.enty[i].sxone >> theCurrentTemp.enty[i].qmin >> theCurrentTemp.enty[i].clslenx;
440 
441  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 2, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
442 
443  for (j=0; j<2; ++j) {
444 
445  db >> theCurrentTemp.enty[i].xpar[j][0] >> theCurrentTemp.enty[i].xpar[j][1]
446  >> theCurrentTemp.enty[i].xpar[j][2] >> theCurrentTemp.enty[i].xpar[j][3] >> theCurrentTemp.enty[i].xpar[j][4];
447 
448  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 6, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
449 
450  }
451 
452  qavg_avg = 0.f;
453  for (j=0; j<9; ++j) {
454 
455  for (k=0; k<TSXSIZE; ++k) {db >> theCurrentTemp.enty[i].xtemp[j][k]; qavg_avg += theCurrentTemp.enty[i].xtemp[j][k];}
456 
457  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 7, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
458  }
459  theCurrentTemp.enty[i].qavg_avg = qavg_avg/9.;
460 
461  for (j=0; j<4; ++j) {
462 
463  db >> theCurrentTemp.enty[i].xavg[j] >> theCurrentTemp.enty[i].xrms[j] >> theCurrentTemp.enty[i].xgx0[j] >> theCurrentTemp.enty[i].xgsig[j];
464 
465  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 10, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
466  }
467 
468  for (j=0; j<4; ++j) {
469 
470  db >> theCurrentTemp.enty[i].xflpar[j][0] >> theCurrentTemp.enty[i].xflpar[j][1] >> theCurrentTemp.enty[i].xflpar[j][2]
471  >> theCurrentTemp.enty[i].xflpar[j][3] >> theCurrentTemp.enty[i].xflpar[j][4] >> theCurrentTemp.enty[i].xflpar[j][5];
472 
473  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 11, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
474  }
475 
476  for (j=0; j<4; ++j) {
477 
478  db >> theCurrentTemp.enty[i].chi2xavg[j] >> theCurrentTemp.enty[i].chi2xmin[j] >> theCurrentTemp.enty[i].chi2xavgc2m[j] >> theCurrentTemp.enty[i].chi2xminc2m[j];
479 
480  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 12, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
481  }
482 
483  for (j=0; j<4; ++j) {
484 
485  db >> theCurrentTemp.enty[i].xavgc2m[j] >> theCurrentTemp.enty[i].xrmsc2m[j] >> theCurrentTemp.enty[i].xgx0c2m[j] >> theCurrentTemp.enty[i].xgsigc2m[j];
486 
487  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
488  }
489 
490  for (j=0; j<4; ++j) {
491 
492  db >> theCurrentTemp.enty[i].xavggen[j] >> theCurrentTemp.enty[i].xrmsgen[j] >> theCurrentTemp.enty[i].xgx0gen[j] >> theCurrentTemp.enty[i].xgsiggen[j];
493 
494  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14b, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
495  }
496 
497  for (j=0; j<4; ++j) {
498 
499  db >> theCurrentTemp.enty[i].xavgbcn[j] >> theCurrentTemp.enty[i].xrmsbcn[j] >> theCurrentTemp.enty[i].xgx0bcn[j] >> theCurrentTemp.enty[i].xgsigbcn[j];
500 
501  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 14c, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
502  }
503 
504  db >> theCurrentTemp.enty[i].chi2xavgone >> theCurrentTemp.enty[i].chi2xminone >> theCurrentTemp.enty[i].qmin2
505  >> theCurrentTemp.enty[i].mpvvav >> theCurrentTemp.enty[i].sigmavav >> theCurrentTemp.enty[i].kappavav
506  >> theCurrentTemp.enty[i].mpvvav2 >> theCurrentTemp.enty[i].sigmavav2 >> theCurrentTemp.enty[i].kappavav2 >> theCurrentTemp.enty[i].spare[0];
507 
508  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 15, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
509 
510  db >> theCurrentTemp.enty[i].qbfrac[0] >> theCurrentTemp.enty[i].qbfrac[1] >> theCurrentTemp.enty[i].qbfrac[2] >> theCurrentTemp.enty[i].fracxone
511  >> theCurrentTemp.enty[i].spare[1] >> theCurrentTemp.enty[i].spare[2] >> theCurrentTemp.enty[i].spare[3] >> theCurrentTemp.enty[i].spare[4]
512  >> theCurrentTemp.enty[i].spare[5] >> theCurrentTemp.enty[i].spare[6];
513 
514  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 16, no template load, run # " << theCurrentTemp.enty[i].runnum << ENDL; return false;}
515 
516  }
517 
518  // next, loop over all barrel x-angle entries
519 
520  for (k=0; k < theCurrentTemp.head.NTyx; ++k) {
521 
522  for (i=0; i < theCurrentTemp.head.NTxx; ++i) {
523 
524  db >> theCurrentTemp.entx[k][i].runnum >> theCurrentTemp.entx[k][i].costrk[0]
525  >> theCurrentTemp.entx[k][i].costrk[1] >> theCurrentTemp.entx[k][i].costrk[2];
526 
527  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 17, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
528 
529  // Calculate the alpha, beta, and cot(beta) for this entry
530 
531  theCurrentTemp.entx[k][i].alpha = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[0]));
532 
533  theCurrentTemp.entx[k][i].cotalpha = theCurrentTemp.entx[k][i].costrk[0]/theCurrentTemp.entx[k][i].costrk[2];
534 
535  theCurrentTemp.entx[k][i].beta = static_cast<float>(atan2((double)theCurrentTemp.entx[k][i].costrk[2], (double)theCurrentTemp.entx[k][i].costrk[1]));
536 
537  theCurrentTemp.entx[k][i].cotbeta = theCurrentTemp.entx[k][i].costrk[1]/theCurrentTemp.entx[k][i].costrk[2];
538 
539  db >> theCurrentTemp.entx[k][i].qavg >> theCurrentTemp.entx[k][i].sxmax >> theCurrentTemp.entx[k][i].dxone >> theCurrentTemp.entx[k][i].sxone >> theCurrentTemp.entx[k][i].qmin >> theCurrentTemp.entx[k][i].clslenx;
540 
541  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 18, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
542 
543  for (j=0; j<2; ++j) {
544 
545  db >> theCurrentTemp.entx[k][i].xpar[j][0] >> theCurrentTemp.entx[k][i].xpar[j][1]
546  >> theCurrentTemp.entx[k][i].xpar[j][2] >> theCurrentTemp.entx[k][i].xpar[j][3] >> theCurrentTemp.entx[k][i].xpar[j][4];
547 
548  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 19, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
549 
550  }
551 
552  qavg_avg = 0.f;
553  for (j=0; j<9; ++j) {
554 
555  for (l=0; l<TSXSIZE; ++k) {db >> theCurrentTemp.entx[k][i].xtemp[j][l]; qavg_avg += theCurrentTemp.entx[k][i].xtemp[j][l];}
556 
557  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 20, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
558  }
559  theCurrentTemp.entx[k][i].qavg_avg = qavg_avg/9.;
560 
561  for (j=0; j<4; ++j) {
562 
563  db >> theCurrentTemp.entx[k][i].xavg[j] >> theCurrentTemp.entx[k][i].xrms[j] >> theCurrentTemp.entx[k][i].xgx0[j] >> theCurrentTemp.entx[k][i].xgsig[j];
564 
565  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 21, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
566  }
567 
568  for (j=0; j<4; ++j) {
569 
570  db >> theCurrentTemp.entx[k][i].xflpar[j][0] >> theCurrentTemp.entx[k][i].xflpar[j][1] >> theCurrentTemp.entx[k][i].xflpar[j][2]
571  >> theCurrentTemp.entx[k][i].xflpar[j][3] >> theCurrentTemp.entx[k][i].xflpar[j][4] >> theCurrentTemp.entx[k][i].xflpar[j][5];
572 
573  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 22, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
574  }
575 
576  for (j=0; j<4; ++j) {
577 
578  db >> theCurrentTemp.entx[k][i].chi2xavg[j] >> theCurrentTemp.entx[k][i].chi2xmin[j] >> theCurrentTemp.entx[k][i].chi2xavgc2m[j] >> theCurrentTemp.entx[k][i].chi2xminc2m[j];
579 
580  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 23, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
581  }
582 
583  for (j=0; j<4; ++j) {
584 
585  db >> theCurrentTemp.entx[k][i].xavgc2m[j] >> theCurrentTemp.entx[k][i].xrmsc2m[j] >> theCurrentTemp.entx[k][i].xgx0c2m[j] >> theCurrentTemp.entx[k][i].xgsigc2m[j];
586 
587  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 24, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
588  }
589 
590  for (j=0; j<4; ++j) {
591 
592  db >> theCurrentTemp.entx[k][i].xavggen[j] >> theCurrentTemp.entx[k][i].xrmsgen[j] >> theCurrentTemp.entx[k][i].xgx0gen[j] >> theCurrentTemp.entx[k][i].xgsiggen[j];
593 
594  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 25, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
595  }
596 
597  for (j=0; j<4; ++j) {
598 
599  db >> theCurrentTemp.entx[k][i].xavgbcn[j] >> theCurrentTemp.entx[k][i].xrmsbcn[j] >> theCurrentTemp.entx[k][i].xgx0bcn[j] >> theCurrentTemp.entx[k][i].xgsigbcn[j];
600 
601  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 26, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
602  }
603 
604  db >> theCurrentTemp.entx[k][i].chi2xavgone >> theCurrentTemp.entx[k][i].chi2xminone >> theCurrentTemp.entx[k][i].qmin2
605  >> theCurrentTemp.entx[k][i].mpvvav >> theCurrentTemp.entx[k][i].sigmavav >> theCurrentTemp.entx[k][i].kappavav
606  >> theCurrentTemp.entx[k][i].mpvvav2 >> theCurrentTemp.entx[k][i].sigmavav2 >> theCurrentTemp.entx[k][i].kappavav2 >> theCurrentTemp.entx[k][i].spare[0];
607 
608  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 27, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
609 
610  db >> theCurrentTemp.entx[k][i].qbfrac[0] >> theCurrentTemp.entx[k][i].qbfrac[1] >> theCurrentTemp.entx[k][i].qbfrac[2] >> theCurrentTemp.entx[k][i].fracxone
611  >> theCurrentTemp.entx[k][i].spare[1] >> theCurrentTemp.entx[k][i].spare[2] >> theCurrentTemp.entx[k][i].spare[3] >> theCurrentTemp.entx[k][i].spare[4]
612  >> theCurrentTemp.entx[k][i].spare[5] >> theCurrentTemp.entx[k][i].spare[6];
613 
614  if(db.fail()) {LOGERROR("SiStripTemplate") << "Error reading file 28, no template load, run # " << theCurrentTemp.entx[k][i].runnum << ENDL; return false;}
615 
616 
617 
618  }
619  }
620 
621 
622  // Add this template to the store
623 
624  theStripTemp_.push_back(theCurrentTemp);
625 
626  }
627  return true;
628 
629 } // TempInit
630 
631 #endif
632 
633 
634 // ************************************************************************************************************
640 // ************************************************************************************************************
641 bool SiStripTemplate::interpolate(int id, float cotalpha, float cotbeta, float locBy)
642 {
643  // Interpolate for a new set of track angles
644 
645  // Local variables
646  int i, j;
647  int ilow, ihigh, iylow, iyhigh, Ny, Nxx, Nyx, imidy, imaxx;
648  float yratio, yxratio, xxratio, sxmax, qcorrect, qxtempcor, chi2xavgone, chi2xminone, cota, cotb, cotalpha0, cotbeta0;
649  bool flip_x;
650 // std::vector <float> xrms(4), xgsig(4), xrmsc2m(4), xgsigc2m(4);
651  std::vector <float> chi2xavg(4), chi2xmin(4), chi2xavgc2m(4), chi2xminc2m(4);
652 
653 
654 // Check to see if interpolation is valid
655 
656 if(id != id_current_ || cotalpha != cota_current_ || cotbeta != cotb_current_) {
657 
658  cota_current_ = cotalpha; cotb_current_ = cotbeta; success_ = true;
659 
660  if(id != id_current_) {
661 
662 // Find the index corresponding to id
663 
664  index_id_ = -1;
665  for(i=0; i<(int)theStripTemp_.size(); ++i) {
666 
667  if(id == theStripTemp_[i].head.ID) {
668 
669  index_id_ = i;
670  id_current_ = id;
671 
672 // Copy the charge scaling factor to the private variable
673 
674  qscale_ = theStripTemp_[index_id_].head.qscale;
675 
676 // Copy the pseudopixel signal size to the private variable
677 
678  s50_ = theStripTemp_[index_id_].head.s50;
679 
680 // Pixel sizes to the private variables
681 
682  xsize_ = theStripTemp_[index_id_].head.xsize;
683  ysize_ = theStripTemp_[index_id_].head.ysize;
684  zsize_ = theStripTemp_[index_id_].head.zsize;
685 
686  break;
687  }
688  }
689  }
690 
691 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
692  if(index_id_ < 0 || index_id_ >= (int)theStripTemp_.size()) {
693  throw cms::Exception("DataCorrupt") << "SiStripTemplate::interpolate can't find needed template ID = " << id << std::endl;
694  }
695 #else
696  assert(index_id_ >= 0 && index_id_ < (int)theStripTemp_.size());
697 #endif
698 
699 // Interpolate the absolute value of cot(beta)
700 
701  abs_cotb_ = std::abs(cotbeta);
702  cotb = abs_cotb_;
703 
704 // qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
705 
706  cotalpha0 = theStripTemp_[index_id_].enty[0].cotalpha;
707  qcorrect=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta*cotbeta+cotalpha0*cotalpha0));
708 // flip quantities when the magnetic field in in the positive y local direction
709  if(locBy > 0.f) {
710  flip_x = true;
711  } else {
712  flip_x = false;
713  }
714 
715  Ny = theStripTemp_[index_id_].head.NTy;
716  Nyx = theStripTemp_[index_id_].head.NTyx;
717  Nxx = theStripTemp_[index_id_].head.NTxx;
718 
719 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
720  if(Ny < 2 || Nyx < 1 || Nxx < 2) {
721  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
722  }
723 #else
724  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
725 #endif
726  imaxx = Nyx - 1;
727  imidy = Nxx/2;
728 
729 // next, loop over all y-angle entries
730 
731  ilow = 0;
732  yratio = 0.f;
733 
734  if(cotb >= theStripTemp_[index_id_].enty[Ny-1].cotbeta) {
735 
736  ilow = Ny-2;
737  yratio = 1.f;
738  success_ = false;
739 
740  } else {
741 
742  if(cotb >= theStripTemp_[index_id_].enty[0].cotbeta) {
743 
744  for (i=0; i<Ny-1; ++i) {
745 
746  if( theStripTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < theStripTemp_[index_id_].enty[i+1].cotbeta) {
747 
748  ilow = i;
749  yratio = (cotb - theStripTemp_[index_id_].enty[i].cotbeta)/(theStripTemp_[index_id_].enty[i+1].cotbeta - theStripTemp_[index_id_].enty[i].cotbeta);
750  break;
751  }
752  }
753  } else { success_ = false; }
754  }
755 
756  ihigh=ilow + 1;
757 
758 // Interpolate/store all y-related quantities (flip displacements when flip_y)
759 
760  yratio_ = yratio;
761  qavg_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].qavg + yratio*theStripTemp_[index_id_].enty[ihigh].qavg;
762  qavg_ *= qcorrect;
763  sxmax = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].sxmax + yratio*theStripTemp_[index_id_].enty[ihigh].sxmax;
764  syparmax_ = sxmax;
765  qmin_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].qmin + yratio*theStripTemp_[index_id_].enty[ihigh].qmin;
766  qmin_ *= qcorrect;
767  qmin2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].qmin2 + yratio*theStripTemp_[index_id_].enty[ihigh].qmin2;
768  qmin2_ *= qcorrect;
769  mpvvav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].mpvvav + yratio*theStripTemp_[index_id_].enty[ihigh].mpvvav;
770  mpvvav_ *= qcorrect;
771  sigmavav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].sigmavav + yratio*theStripTemp_[index_id_].enty[ihigh].sigmavav;
772  kappavav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].kappavav + yratio*theStripTemp_[index_id_].enty[ihigh].kappavav;
773  mpvvav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].mpvvav2 + yratio*theStripTemp_[index_id_].enty[ihigh].mpvvav2;
774  mpvvav2_ *= qcorrect;
775  sigmavav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].sigmavav2 + yratio*theStripTemp_[index_id_].enty[ihigh].sigmavav2;
776  kappavav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].kappavav2 + yratio*theStripTemp_[index_id_].enty[ihigh].kappavav2;
777  qavg_avg_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].qavg_avg + yratio*theStripTemp_[index_id_].enty[ihigh].qavg_avg;
778  qavg_avg_ *= qcorrect;
779  for(i=0; i<2 ; ++i) {
780  for(j=0; j<5 ; ++j) {
781 // Charge loss switches sides when cot(alpha) changes sign
782  if(flip_x) {
783  xparly0_[1-i][j] = theStripTemp_[index_id_].enty[ilow].xpar[i][j];
784  xparhy0_[1-i][j] = theStripTemp_[index_id_].enty[ihigh].xpar[i][j];
785  } else {
786  xparly0_[i][j] = theStripTemp_[index_id_].enty[ilow].xpar[i][j];
787  xparhy0_[i][j] = theStripTemp_[index_id_].enty[ihigh].xpar[i][j];
788  }
789  }
790  }
791  for(i=0; i<4; ++i) {
792  chi2xavg[i]=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xavg[i] + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xavg[i];
793  chi2xmin[i]=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xmin[i] + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xmin[i];
794  chi2xavgc2m[i]=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xavgc2m[i] + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xavgc2m[i];
795  chi2xminc2m[i]=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xminc2m[i] + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xminc2m[i];
796  }
797 
799 
800  chi2xavgone=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xavgone + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xavgone;
801  chi2xminone=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].chi2xminone + yratio*theStripTemp_[index_id_].enty[ihigh].chi2xminone;
802  // for(i=0; i<10; ++i) {
803 // pyspare[i]=(1.f - yratio)*theStripTemp_[index_id_].enty[ilow].yspare[i] + yratio*theStripTemp_[index_id_].enty[ihigh].yspare[i];
804 // }
805 
806 
807 // next, loop over all x-angle entries, first, find relevant y-slices
808 
809  iylow = 0;
810  yxratio = 0.f;
811 
812  if(abs_cotb_ >= theStripTemp_[index_id_].entx[Nyx-1][0].cotbeta) {
813 
814  iylow = Nyx-2;
815  yxratio = 1.f;
816 
817  } else if(abs_cotb_ >= theStripTemp_[index_id_].entx[0][0].cotbeta) {
818 
819  for (i=0; i<Nyx-1; ++i) {
820 
821  if( theStripTemp_[index_id_].entx[i][0].cotbeta <= abs_cotb_ && abs_cotb_ < theStripTemp_[index_id_].entx[i+1][0].cotbeta) {
822 
823  iylow = i;
824  yxratio = (abs_cotb_ - theStripTemp_[index_id_].entx[i][0].cotbeta)/(theStripTemp_[index_id_].entx[i+1][0].cotbeta - theStripTemp_[index_id_].entx[i][0].cotbeta);
825  break;
826  }
827  }
828  }
829 
830  iyhigh=iylow + 1;
831 
832  ilow = 0;
833  xxratio = 0.f;
834  if(flip_x) {cota = -cotalpha;} else {cota = cotalpha;}
835 
836  if(cota >= theStripTemp_[index_id_].entx[0][Nxx-1].cotalpha) {
837 
838  ilow = Nxx-2;
839  xxratio = 1.f;
840  success_ = false;
841 
842  } else {
843 
844  if(cota >= theStripTemp_[index_id_].entx[0][0].cotalpha) {
845 
846  for (i=0; i<Nxx-1; ++i) {
847 
848  if( theStripTemp_[index_id_].entx[0][i].cotalpha <= cota && cota < theStripTemp_[index_id_].entx[0][i+1].cotalpha) {
849 
850  ilow = i;
851  xxratio = (cota - theStripTemp_[index_id_].entx[0][i].cotalpha)/(theStripTemp_[index_id_].entx[0][i+1].cotalpha - theStripTemp_[index_id_].entx[0][i].cotalpha);
852  break;
853  }
854  }
855  } else { success_ = false; }
856  }
857 
858  ihigh=ilow + 1;
859 
860 // Interpolate/store all x-related quantities
861 
862  yxratio_ = yxratio;
863  xxratio_ = xxratio;
864 
865 // sxparmax defines the maximum charge for which the parameters xpar are defined (not rescaled by cotbeta)
866 
867  sxparmax_ = (1.f - xxratio)*theStripTemp_[index_id_].entx[imaxx][ilow].sxmax + xxratio*theStripTemp_[index_id_].entx[imaxx][ihigh].sxmax;
868  sxmax_ = sxparmax_;
869  if(theStripTemp_[index_id_].entx[imaxx][imidy].sxmax != 0.f) {sxmax_=sxmax_/theStripTemp_[index_id_].entx[imaxx][imidy].sxmax*sxmax;}
870  dxone_ = (1.f - xxratio)*theStripTemp_[index_id_].entx[0][ilow].dxone + xxratio*theStripTemp_[index_id_].entx[0][ihigh].dxone;
871  if(flip_x) {dxone_ = -dxone_;}
872  sxone_ = (1.f - xxratio)*theStripTemp_[index_id_].entx[0][ilow].sxone + xxratio*theStripTemp_[index_id_].entx[0][ihigh].sxone;
873  clslenx_ = fminf(theStripTemp_[index_id_].entx[0][ilow].clslenx, theStripTemp_[index_id_].entx[0][ihigh].clslenx);
874  for(i=0; i<2 ; ++i) {
875  for(j=0; j<5 ; ++j) {
876  if(flip_x) {
877  xpar0_[1-i][j] = theStripTemp_[index_id_].entx[imaxx][imidy].xpar[i][j];
878  xparl_[1-i][j] = theStripTemp_[index_id_].entx[imaxx][ilow].xpar[i][j];
879  xparh_[1-i][j] = theStripTemp_[index_id_].entx[imaxx][ihigh].xpar[i][j];
880  } else {
881  xpar0_[i][j] = theStripTemp_[index_id_].entx[imaxx][imidy].xpar[i][j];
882  xparl_[i][j] = theStripTemp_[index_id_].entx[imaxx][ilow].xpar[i][j];
883  xparh_[i][j] = theStripTemp_[index_id_].entx[imaxx][ihigh].xpar[i][j];
884  }
885  }
886  }
887 
888 // sxmax is the maximum allowed strip charge (used for truncation)
889 
890  sxmax_=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].sxmax + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].sxmax)
891  +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].sxmax + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].sxmax);
892 
893  for(i=0; i<4; ++i) {
894  xavg_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xavg[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xavg[i])
895  +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xavg[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xavg[i]);
896  if(flip_x) {xavg_[i] = -xavg_[i];}
897 
898  xrms_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xrms[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xrms[i])
899  +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xrms[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xrms[i]);
900 
901 // xgx0_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xgx0[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xgx0[i])
902 // +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xgx0[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xgx0[i]);
903 
904 // xgsig_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xgsig[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xgsig[i])
905 // +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xgsig[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xgsig[i]);
906 
907  xavgc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xavgc2m[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xavgc2m[i])
908  +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xavgc2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xavgc2m[i]);
909  if(flip_x) {xavgc2m_[i] = -xavgc2m_[i];}
910 
911  xrmsc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xrmsc2m[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xrmsc2m[i])
912  +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xrmsc2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xrmsc2m[i]);
913 
914  xavgbcn_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xavgbcn[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xavgbcn[i])
915  +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xavgbcn[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xavgbcn[i]);
916  if(flip_x) {xavgbcn_[i] = -xavgbcn_[i];}
917 
918  xrmsbcn_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xrmsbcn[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xrmsbcn[i])
919  +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xrmsbcn[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xrmsbcn[i]);
920 
921 // xgx0c2m_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xgx0c2m[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xgx0c2m[i])
922 // +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xgx0c2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xgx0c2m[i]);
923 
924 // xgsigc2m_[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xgsigc2m[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xgsigc2m[i])
925 // +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xgsigc2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xgsigc2m[i]);
926 //
927 // Try new interpolation scheme
928 //
929 // chi2xavg_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[imaxx][ilow].chi2xavg[i] + xxratio*theStripTemp_[index_id_].entx[imaxx][ihigh].chi2xavg[i]);
930 // if(theStripTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i] != 0.f) {chi2xavg_[i]=chi2xavg_[i]/theStripTemp_[index_id_].entx[imaxx][imidy].chi2xavg[i]*chi2xavg[i];}
931 //
932 // chi2xmin_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[imaxx][ilow].chi2xmin[i] + xxratio*theStripTemp_[index_id_].entx[imaxx][ihigh].chi2xmin[i]);
933 // if(theStripTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i] != 0.f) {chi2xmin_[i]=chi2xmin_[i]/theStripTemp_[index_id_].entx[imaxx][imidy].chi2xmin[i]*chi2xmin[i];}
934 //
935  chi2xavg_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xavg[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xavg[i]);
936  if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavg[i] != 0.f) {chi2xavg_[i]=chi2xavg_[i]/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavg[i]*chi2xavg[i];}
937 
938  chi2xmin_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xmin[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xmin[i]);
939  if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xmin[i] != 0.f) {chi2xmin_[i]=chi2xmin_[i]/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xmin[i]*chi2xmin[i];}
940 
941  chi2xavgc2m_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xavgc2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xavgc2m[i]);
942  if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavgc2m[i] != 0.f) {chi2xavgc2m_[i]=chi2xavgc2m_[i]/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavgc2m[i]*chi2xavgc2m[i];}
943 
944  chi2xminc2m_[i]=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xminc2m[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xminc2m[i]);
945  if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xminc2m[i] != 0.f) {chi2xminc2m_[i]=chi2xminc2m_[i]/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xminc2m[i]*chi2xminc2m[i];}
946 
947  for(j=0; j<6 ; ++j) {
948  xflparll_[i][j] = theStripTemp_[index_id_].entx[iylow][ilow].xflpar[i][j];
949  xflparlh_[i][j] = theStripTemp_[index_id_].entx[iylow][ihigh].xflpar[i][j];
950  xflparhl_[i][j] = theStripTemp_[index_id_].entx[iyhigh][ilow].xflpar[i][j];
951  xflparhh_[i][j] = theStripTemp_[index_id_].entx[iyhigh][ihigh].xflpar[i][j];
952  // Since Q_fl is odd under cotbeta, it flips qutomatically, change only even terms
953 
954  if(flip_x && (j == 0 || j == 2 || j == 4)) {
955  xflparll_[i][j] = -xflparll_[i][j];
956  xflparlh_[i][j] = -xflparlh_[i][j];
957  xflparhl_[i][j] = -xflparhl_[i][j];
958  xflparhh_[i][j] = -xflparhh_[i][j];
959  }
960  }
961  }
962 
963 // Do the spares next
964 
965  chi2xavgone_=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xavgone + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xavgone);
966  if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavgone != 0.f) {chi2xavgone_=chi2xavgone_/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xavgone*chi2xavgone;}
967 
968  chi2xminone_=((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].chi2xminone + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].chi2xminone);
969  if(theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xminone != 0.f) {chi2xminone_=chi2xminone_/theStripTemp_[index_id_].entx[iyhigh][imidy].chi2xminone*chi2xminone;}
970  // for(i=0; i<10; ++i) {
971 // pxspare[i]=(1.f - yxratio)*((1.f - xxratio)*theStripTemp_[index_id_].entx[iylow][ilow].xspare[i] + xxratio*theStripTemp_[index_id_].entx[iylow][ihigh].xspare[i])
972 // +yxratio*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xspare[i] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xspare[i]);
973 // }
974 
975 // Interpolate and build the x-template
976 
977 // qxtempcor corrects the total charge to the actual track angles (not actually needed for the template fits, but useful for Guofan)
978 
979  cotbeta0 = theStripTemp_[index_id_].entx[iyhigh][0].cotbeta;
980  qxtempcor=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta0*cotbeta0+cotalpha*cotalpha));
981 
982  for(i=0; i<9; ++i) {
983  xtemp_[i][0] = 0.f;
984  xtemp_[i][1] = 0.f;
985  xtemp_[i][BSXM2] = 0.f;
986  xtemp_[i][BSXM1] = 0.f;
987  for(j=0; j<TSXSIZE; ++j) {
988 // Take next largest x-slice for the x-template (it reduces bias in the forward direction after irradiation)
989 // xtemp_[i][j+2]=(1.f - xxratio)*theStripTemp_[index_id_].entx[imaxx][ilow].xtemp[i][j] + xxratio*theStripTemp_[index_id_].entx[imaxx][ihigh].xtemp[i][j];
990  if(flip_x) {
991  xtemp_[8-i][BSXM3-j]=qxtempcor*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xtemp[i][j] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xtemp[i][j]);
992  } else {
993  xtemp_[i][j+2]=qxtempcor*((1.f - xxratio)*theStripTemp_[index_id_].entx[iyhigh][ilow].xtemp[i][j] + xxratio*theStripTemp_[index_id_].entx[iyhigh][ihigh].xtemp[i][j]);
994  }
995  }
996  }
997 
998  lorxwidth_ = theStripTemp_[index_id_].head.lorxwidth;
999  if(locBy > 0.f) {lorxwidth_ = -lorxwidth_;}
1000 
1001  }
1002 
1003  return success_;
1004 } // interpolate
1005 
1006 
1007 
1008 // ************************************************************************************************************
1016 // ************************************************************************************************************
1017  void SiStripTemplate::xsigma2(int fxpix, int lxpix, float sxthr, float xsum[11], float xsig2[11])
1018 
1019 {
1020  // Interpolate using quantities already stored in the private variables
1021 
1022  // Local variables
1023  int i;
1024  float sigi, sigi2, sigi3, sigi4, yint, sxmax, x0, qscale;
1025  float sigiy, sigiy2, sigiy3, sigiy4;
1026 
1027  // Make sure that input is OK
1028 
1029 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1030  if(fxpix < 2 || fxpix >= BSXM2) {
1031  throw cms::Exception("DataCorrupt") << "SiStripTemplate::xsigma2 called with fxpix = " << fxpix << std::endl;
1032  }
1033 #else
1034  assert(fxpix > 1 && fxpix < BSXM2);
1035 #endif
1036 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1037  if(lxpix < fxpix || lxpix >= BSXM2) {
1038  throw cms::Exception("DataCorrupt") << "SiStripTemplate::xsigma2 called with lxpix/fxpix = " << lxpix << "/" << fxpix << std::endl;
1039  }
1040 #else
1041  assert(lxpix >= fxpix && lxpix < BSXM2);
1042 #endif
1043 
1044 // Define the maximum signal to use in the parameterization
1045 
1046  sxmax = sxmax_;
1047  if(sxmax_ > sxparmax_) {sxmax = sxparmax_;}
1048 
1049 // Evaluate pixel-by-pixel uncertainties (weights) for the templ analysis
1050 
1051  for(i=fxpix-2; i<=lxpix+2; ++i) {
1052  if(i < fxpix || i > lxpix) {
1053 
1054 // Nearest pseudopixels have uncertainties of 50% of threshold, next-nearest have 10% of threshold
1055 
1056  xsig2[i] = s50_*s50_;
1057  } else {
1058  if(xsum[i] < sxmax) {
1059  sigi = xsum[i];
1060  qscale = 1.f;
1061  } else {
1062  sigi = sxmax;
1063  qscale = xsum[i]/sxmax;
1064  }
1065  sigi2 = sigi*sigi; sigi3 = sigi2*sigi; sigi4 = sigi3*sigi;
1066 
1067  if(xsum[i] < syparmax_) {
1068  sigiy = xsum[i];
1069  } else {
1070  sigiy = syparmax_;
1071  }
1072  sigiy2 = sigiy*sigiy; sigiy3 = sigiy2*sigiy; sigiy4 = sigiy3*sigiy;
1073 
1074 // First, do the cotbeta interpolation
1075 
1076  if(i <= BSHX) {
1077  yint = (1.f-yratio_)*
1078  (xparly0_[0][0]+xparly0_[0][1]*sigiy+xparly0_[0][2]*sigiy2+xparly0_[0][3]*sigiy3+xparly0_[0][4]*sigiy4)
1079  + yratio_*
1080  (xparhy0_[0][0]+xparhy0_[0][1]*sigiy+xparhy0_[0][2]*sigiy2+xparhy0_[0][3]*sigiy3+xparhy0_[0][4]*sigiy4);
1081  } else {
1082  yint = (1.f-yratio_)*
1083  (xparly0_[1][0]+xparly0_[1][1]*sigiy+xparly0_[1][2]*sigiy2+xparly0_[1][3]*sigiy3+xparly0_[1][4]*sigiy4)
1084  + yratio_*
1085  (xparhy0_[1][0]+xparhy0_[1][1]*sigiy+xparhy0_[1][2]*sigiy2+xparhy0_[1][3]*sigiy3+xparhy0_[1][4]*sigiy4);
1086  }
1087 
1088 // Next, do the cotalpha interpolation
1089 
1090  if(i <= BSHX) {
1091  xsig2[i] = (1.f-xxratio_)*
1092  (xparl_[0][0]+xparl_[0][1]*sigi+xparl_[0][2]*sigi2+xparl_[0][3]*sigi3+xparl_[0][4]*sigi4)
1093  + xxratio_*
1094  (xparh_[0][0]+xparh_[0][1]*sigi+xparh_[0][2]*sigi2+xparh_[0][3]*sigi3+xparh_[0][4]*sigi4);
1095  } else {
1096  xsig2[i] = (1.f-xxratio_)*
1097  (xparl_[1][0]+xparl_[1][1]*sigi+xparl_[1][2]*sigi2+xparl_[1][3]*sigi3+xparl_[1][4]*sigi4)
1098  + xxratio_*
1099  (xparh_[1][0]+xparh_[1][1]*sigi+xparh_[1][2]*sigi2+xparh_[1][3]*sigi3+xparh_[1][4]*sigi4);
1100  }
1101 
1102 // Finally, get the mid-point value of the cotalpha function
1103 
1104  if(i <= BSHX) {
1105  x0 = xpar0_[0][0]+xpar0_[0][1]*sigi+xpar0_[0][2]*sigi2+xpar0_[0][3]*sigi3+xpar0_[0][4]*sigi4;
1106  } else {
1107  x0 = xpar0_[1][0]+xpar0_[1][1]*sigi+xpar0_[1][2]*sigi2+xpar0_[1][3]*sigi3+xpar0_[1][4]*sigi4;
1108  }
1109 
1110 // Finally, rescale the yint value for cotalpha variation
1111 
1112  if(x0 != 0.f) {xsig2[i] = xsig2[i]/x0 * yint;}
1113  xsig2[i] *=qscale;
1114  if(xsum[i] > sxthr) {xsig2[i] = 1.e8f;}
1115  if(xsig2[i] <= 0.f) {LOGERROR("SiStripTemplate") << "neg x-error-squared = " << xsig2[i] << ", id = " << id_current_ << ", index = " << index_id_ <<
1116  ", cot(alpha) = " << cota_current_ << ", cot(beta) = " << cotb_current_ << ", sigi = " << sigi << ", sxparmax = " << sxparmax_ << ", sxmax = " << sxmax_ << ENDL;}
1117  }
1118  }
1119 
1120  return;
1121 
1122 } // End xsigma2
1123 
1124 
1125 
1126 
1127 // ************************************************************************************************************
1131 // ************************************************************************************************************
1132  float SiStripTemplate::xflcorr(int binq, float qflx)
1133 
1134 {
1135  // Interpolate using quantities already stored in the private variables
1136 
1137  // Local variables
1138  float qfl, qfl2, qfl3, qfl4, qfl5, dx;
1139 
1140  // Make sure that input is OK
1141 
1142 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1143  if(binq < 0 || binq > 3) {
1144  throw cms::Exception("DataCorrupt") << "SiStripTemplate::xflcorr called with binq = " << binq << std::endl;
1145  }
1146 #else
1147  assert(binq >= 0 && binq < 4);
1148 #endif
1149 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1150  if(fabs((double)qflx) > 1.) {
1151  throw cms::Exception("DataCorrupt") << "SiStripTemplate::xflcorr called with qflx = " << qflx << std::endl;
1152  }
1153 #else
1154  assert(fabs((double)qflx) <= 1.);
1155 #endif
1156 
1157 // Define the maximum signal to allow before de-weighting a pixel
1158 
1159  qfl = qflx;
1160 
1161  if(qfl < -0.9f) {qfl = -0.9f;}
1162  if(qfl > 0.9f) {qfl = 0.9f;}
1163 
1164 // Interpolate between the two polynomials
1165 
1166  qfl2 = qfl*qfl; qfl3 = qfl2*qfl; qfl4 = qfl3*qfl; qfl5 = qfl4*qfl;
1167  dx = (1.f - yxratio_)*((1.f-xxratio_)*(xflparll_[binq][0]+xflparll_[binq][1]*qfl+xflparll_[binq][2]*qfl2+xflparll_[binq][3]*qfl3+xflparll_[binq][4]*qfl4+xflparll_[binq][5]*qfl5)
1168  + xxratio_*(xflparlh_[binq][0]+xflparlh_[binq][1]*qfl+xflparlh_[binq][2]*qfl2+xflparlh_[binq][3]*qfl3+xflparlh_[binq][4]*qfl4+xflparlh_[binq][5]*qfl5))
1169  + yxratio_*((1.f-xxratio_)*(xflparhl_[binq][0]+xflparhl_[binq][1]*qfl+xflparhl_[binq][2]*qfl2+xflparhl_[binq][3]*qfl3+xflparhl_[binq][4]*qfl4+xflparhl_[binq][5]*qfl5)
1170  + xxratio_*(xflparhh_[binq][0]+xflparhh_[binq][1]*qfl+xflparhh_[binq][2]*qfl2+xflparhh_[binq][3]*qfl3+xflparhh_[binq][4]*qfl4+xflparhh_[binq][5]*qfl5));
1171 
1172  return dx;
1173 
1174 } // End xflcorr
1175 
1176 
1177 // ************************************************************************************************************
1182 // ************************************************************************************************************
1183  void SiStripTemplate::xtemp(int fxbin, int lxbin, float xtemplate[41][BSXSIZE])
1184 
1185 {
1186  // Retrieve already interpolated quantities
1187 
1188  // Local variables
1189  int i, j;
1190 
1191  // Verify that input parameters are in valid range
1192 
1193 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1194  if(fxbin < 0 || fxbin > 40) {
1195  throw cms::Exception("DataCorrupt") << "SiStripTemplate::xtemp called with fxbin = " << fxbin << std::endl;
1196  }
1197 #else
1198  assert(fxbin >= 0 && fxbin < 41);
1199 #endif
1200 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1201  if(lxbin < 0 || lxbin > 40) {
1202  throw cms::Exception("DataCorrupt") << "SiStripTemplate::xtemp called with lxbin = " << lxbin << std::endl;
1203  }
1204 #else
1205  assert(lxbin >= 0 && lxbin < 41);
1206 #endif
1207 
1208 // Build the x-template, the central 25 bins are here in all cases
1209 
1210  for(i=0; i<9; ++i) {
1211  for(j=0; j<BSXSIZE; ++j) {
1212  xtemplate[i+16][j]=xtemp_[i][j];
1213  }
1214  }
1215  for(i=0; i<8; ++i) {
1216  xtemplate[i+8][BSXM1] = 0.f;
1217  for(j=0; j<BSXM1; ++j) {
1218  xtemplate[i+8][j]=xtemp_[i][j+1];
1219  }
1220  }
1221  for(i=1; i<9; ++i) {
1222  xtemplate[i+24][0] = 0.f;
1223  for(j=0; j<BSXM1; ++j) {
1224  xtemplate[i+24][j+1]=xtemp_[i][j];
1225  }
1226  }
1227 
1228 // Add more bins if needed
1229 
1230  if(fxbin < 8) {
1231  for(i=0; i<8; ++i) {
1232  xtemplate[i][BSXM2] = 0.f;
1233  xtemplate[i][BSXM1] = 0.f;
1234  for(j=0; j<BSXM2; ++j) {
1235  xtemplate[i][j]=xtemp_[i][j+2];
1236  }
1237  }
1238  }
1239  if(lxbin > 32) {
1240  for(i=1; i<9; ++i) {
1241  xtemplate[i+32][0] = 0.f;
1242  xtemplate[i+32][1] = 0.f;
1243  for(j=0; j<BSXM2; ++j) {
1244  xtemplate[i+32][j+2]=xtemp_[i][j];
1245  }
1246  }
1247  }
1248 
1249  return;
1250 
1251 } // End xtemp
1252 
1253 
1254 // ************************************************************************************************************
1258 // ************************************************************************************************************
1259 void SiStripTemplate::sxtemp(float xhit, std::vector<float>& cluster)
1260 
1261 {
1262  // Retrieve already interpolated quantities
1263 
1264  // Local variables
1265  int i, j;
1266 
1267  // Extract x template based upon the hit position
1268 
1269  float xpix = xhit/xsize_ + 0.5f;
1270 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1271  if(xpix < 0.f) {
1272  throw cms::Exception("DataCorrupt") << "SiStripTemplate::2xtemp called with xhit = " << xhit << std::endl;
1273  }
1274 #else
1275  assert(xpix >= 0.f);
1276 #endif
1277 
1278 // cpix is struck pixel(strip) of the cluster
1279 
1280  int cpix = (int)xpix;
1281  int shift = BSHX - cpix;
1282 
1283 // xbin the floating bin number and cbin is the bin number (between 0 and 7) of the interpolated template
1284 
1285  float xbin = 8.*(xpix-(float)cpix);
1286  int cbin = (int)xbin;
1287 
1288  float xfrac = xbin-(float)cbin;
1289 
1290  int sizex = std::min((int)cluster.size(), BSXSIZE);
1291 
1292 // shift and interpolate the correct cluster shape
1293 
1294  for(i=0; i<sizex; ++i) {
1295  j = i+shift;
1296  if(j < 0 || j > sizex-1) {cluster[i] = 0.f;} else {
1297  cluster[i]=(1.f-xfrac)*xtemp_[cbin][j]+xfrac*xtemp_[cbin+1][j];
1298  if(cluster[i] < s50_) cluster[i] = 0.f;
1299  }
1300 
1301 // Return cluster in same charge units
1302 
1303  cluster[i] /= qscale_;
1304  }
1305 
1306 
1307  return;
1308 
1309 } // End sxtemp
1310 
1311 // ************************************************************************************************************
1313 // ************************************************************************************************************
1314 int SiStripTemplate::cxtemp()
1315 
1316 {
1317  // Retrieve already interpolated quantities
1318 
1319  // Local variables
1320  int j;
1321 
1322  // Analyze only pixels along the central entry
1323  // First, find the maximum signal and then work out to the edges
1324 
1325  float sigmax = 0.f;
1326  int jmax = -1;
1327 
1328  for(j=0; j<BSXSIZE; ++j) {
1329  if(xtemp_[4][j] > sigmax) {
1330  sigmax = xtemp_[4][j];
1331  jmax = j;
1332  }
1333  }
1334  if(sigmax < 2.*s50_ || jmax<1 || jmax>BSXM2) {return -1;}
1335 
1336  // Now search forward and backward
1337 
1338  int jend = jmax;
1339 
1340  for(j=jmax+1; j<BSXM1; ++j) {
1341  if(xtemp_[4][j] < 2.*s50_) break;
1342  jend = j;
1343  }
1344 
1345  int jbeg = jmax;
1346 
1347  for(j=jmax-1; j>0; --j) {
1348  if(xtemp_[4][j] < 2.*s50_) break;
1349  jbeg = j;
1350  }
1351 
1352  return (jbeg+jend)/2;
1353 
1354 } // End cxtemp
1355 
1356 
1357 // ************************************************************************************************************
1361 // ************************************************************************************************************
1362 void SiStripTemplate::xtemp3d_int(int nxpix, int& nxbins)
1363 
1364 {
1365  // Retrieve already interpolated quantities
1366 
1367  // Local variables
1368  int i, j, k;
1369  int ioff0, ioffp, ioffm;
1370 
1371  // Verify that input parameters are in valid range
1372 
1373 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1374  if(nxpix < 1 || nxpix >= BSXM3) {
1375  throw cms::Exception("DataCorrupt") << "SiPixelTemplate::xtemp3d called with nxpix = " << nxpix << std::endl;
1376  }
1377 #else
1378  assert(nxpix > 0 && nxpix < BSXM3);
1379 #endif
1380 
1381  // Calculate the size of the shift in pixels needed to span the entire cluster
1382 
1383  float diff = fabsf(nxpix - clslenx_)/2. + 1.f;
1384  int nshift = (int)diff;
1385  if((diff - nshift) > 0.5f) {++nshift;}
1386 
1387  // Calculate the number of bins needed to specify each hit range
1388 
1389  nxbins_ = 9 + 16*nshift;
1390 
1391  // Create a 2-d working template with the correct size
1392 
1393  temp2dx_.resize(boost::extents[nxbins_][BSXSIZE]);
1394 
1395  // The 9 central bins are copied from the interpolated private store
1396 
1397  ioff0 = 8*nshift;
1398 
1399  for(i=0; i<9; ++i) {
1400  for(j=0; j<BSXSIZE; ++j) {
1401  temp2dx_[i+ioff0][j]=xtemp_[i][j];
1402  }
1403  }
1404 
1405  // Add the +- shifted templates
1406 
1407  for(k=1; k<=nshift; ++k) {
1408  ioffm=ioff0-k*8;
1409  for(i=0; i<8; ++i) {
1410  for(j=0; j<k; ++j) {
1411  temp2dx_[i+ioffm][BSXM1-j] = 0.f;
1412  }
1413  for(j=0; j<BSXSIZE-k; ++j) {
1414  temp2dx_[i+ioffm][j]=xtemp_[i][j+k];
1415  }
1416  }
1417  ioffp=ioff0+k*8;
1418  for(i=1; i<9; ++i) {
1419  for(j=0; j<k; ++j) {
1420  temp2dx_[i+ioffp][j] = 0.f;
1421  }
1422  for(j=0; j<BSXSIZE-k; ++j) {
1423  temp2dx_[i+ioffp][j+k]=xtemp_[i][j];
1424  }
1425  }
1426  }
1427 
1428  nxbins = nxbins_;
1429 
1430  return;
1431 
1432 } // End xtemp3d_int
1433 
1434 
1435 
1436 // ************************************************************************************************************
1440 // ************************************************************************************************************
1441 void SiStripTemplate::xtemp3d(int i, int j, std::vector<float>& xtemplate)
1442 
1443 {
1444  // Sum two 2-d templates to make the 3-d template
1445  if(i >= 0 && i < nxbins_ && j <= i) {
1446  for(int k=0; k<BSXSIZE; ++k) {
1447  xtemplate[k]=temp2dx_[i][k]+temp2dx_[j][k];
1448  }
1449  } else {
1450  for(int k=0; k<BSXSIZE; ++k) {
1451  xtemplate[k]=0.;
1452  }
1453  }
1454 
1455  return;
1456 
1457 } // End xtemp3d
1458 
1459 
1460 
1461 // ************************************************************************************************************
1468 // ************************************************************************************************************
1469 int SiStripTemplate::qbin(int id, float cotalpha, float cotbeta, float qclus)
1470 
1471 {
1472  // Interpolate for a new set of track angles
1473 
1474  // Local variables
1475  int i, binq;
1476  int ilow, ihigh, Ny, Nxx, Nyx, index;
1477  float yratio;
1478  float acotb, qscale, qavg, qmin, qmin2, fq, qtotal, qcorrect, cotalpha0;
1479 
1480 
1481 // Find the index corresponding to id
1482 
1483  index = -1;
1484  for(i=0; i<(int)theStripTemp_.size(); ++i) {
1485 
1486  if(id == theStripTemp_[i].head.ID) {
1487 
1488  index = i;
1489  break;
1490  }
1491  }
1492 
1493 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1494  if(index < 0 || index >= (int)theStripTemp_.size()) {
1495  throw cms::Exception("DataCorrupt") << "SiStripTemplate::qbin can't find needed template ID = " << id << std::endl;
1496  }
1497 #else
1498  assert(index >= 0 && index < (int)theStripTemp_.size());
1499 #endif
1500 
1501 //
1502 
1503 // Interpolate the absolute value of cot(beta)
1504 
1505  acotb = fabs((double)cotbeta);
1506 
1507 // qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
1508 
1509  // qcorrect corrects the cot(alpha)=0 cluster charge for non-zero cot(alpha)
1510 
1511  cotalpha0 = theStripTemp_[index].enty[0].cotalpha;
1512  qcorrect=std::sqrt((1.f+cotbeta*cotbeta+cotalpha*cotalpha)/(1.f+cotbeta*cotbeta+cotalpha0*cotalpha0));
1513 
1514  // Copy the charge scaling factor to the private variable
1515 
1516  qscale = theStripTemp_[index].head.qscale;
1517 
1518  Ny = theStripTemp_[index].head.NTy;
1519  Nyx = theStripTemp_[index].head.NTyx;
1520  Nxx = theStripTemp_[index].head.NTxx;
1521 
1522 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1523  if(Ny < 2 || Nyx < 1 || Nxx < 2) {
1524  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny/Nyx/Nxx = " << Ny << "/" << Nyx << "/" << Nxx << std::endl;
1525  }
1526 #else
1527  assert(Ny > 1 && Nyx > 0 && Nxx > 1);
1528 #endif
1529 
1530 // next, loop over all y-angle entries
1531 
1532  ilow = 0;
1533  yratio = 0.f;
1534 
1535  if(acotb >= theStripTemp_[index].enty[Ny-1].cotbeta) {
1536 
1537  ilow = Ny-2;
1538  yratio = 1.f;
1539 
1540  } else {
1541 
1542  if(acotb >= theStripTemp_[index].enty[0].cotbeta) {
1543 
1544  for (i=0; i<Ny-1; ++i) {
1545 
1546  if( theStripTemp_[index].enty[i].cotbeta <= acotb && acotb < theStripTemp_[index].enty[i+1].cotbeta) {
1547 
1548  ilow = i;
1549  yratio = (acotb - theStripTemp_[index].enty[i].cotbeta)/(theStripTemp_[index].enty[i+1].cotbeta - theStripTemp_[index].enty[i].cotbeta);
1550  break;
1551  }
1552  }
1553  }
1554  }
1555 
1556  ihigh=ilow + 1;
1557 
1558 // Interpolate/store all y-related quantities (flip displacements when flip_y)
1559 
1560  qavg = (1.f - yratio)*theStripTemp_[index].enty[ilow].qavg + yratio*theStripTemp_[index].enty[ihigh].qavg;
1561  qavg *= qcorrect;
1562  qmin = (1.f - yratio)*theStripTemp_[index].enty[ilow].qmin + yratio*theStripTemp_[index].enty[ihigh].qmin;
1563  qmin *= qcorrect;
1564  qmin2 = (1.f - yratio)*theStripTemp_[index].enty[ilow].qmin2 + yratio*theStripTemp_[index].enty[ihigh].qmin2;
1565  qmin2 *= qcorrect;
1566 
1567 
1568 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1569  if(qavg <= 0.f || qmin <= 0.f) {
1570  throw cms::Exception("DataCorrupt") << "SiStripTemplate::qbin, qavg or qmin <= 0,"
1571  << " Probably someone called the generic pixel reconstruction with an illegal trajectory state" << std::endl;
1572  }
1573 #else
1574  assert(qavg > 0.f && qmin > 0.f);
1575 #endif
1576 
1577 // Scale the input charge to account for differences between pixelav and CMSSW simulation or data
1578 
1579  qtotal = qscale*qclus;
1580 
1581 // uncertainty and final corrections depend upon total charge bin
1582 
1583  fq = qtotal/qavg;
1584  if(fq > 1.5f) {
1585  binq=0;
1586  } else {
1587  if(fq > 1.0f) {
1588  binq=1;
1589  } else {
1590  if(fq > 0.85f) {
1591  binq=2;
1592  } else {
1593  binq=3;
1594  }
1595  }
1596  }
1597 
1598 // If the charge is too small (then flag it)
1599 
1600  if(qtotal < 0.95f*qmin) {binq = 5;} else {if(qtotal < 0.95f*qmin2) {binq = 4;}}
1601 
1602  return binq;
1603 
1604 } // qbin
1605 
1606 
1607 // ************************************************************************************************************
1612 // ************************************************************************************************************
1613 void SiStripTemplate::vavilov_pars(double& mpv, double& sigma, double& kappa)
1614 
1615 {
1616  // Local variables
1617  int i;
1618  int ilow, ihigh, Ny;
1619  float yratio, cotb, cotalpha0, arg;
1620 
1621 // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta)
1622 
1623  cotalpha0 = theStripTemp_[index_id_].enty[0].cotalpha;
1624  arg = cotb_current_*cotb_current_ + cota_current_*cota_current_ - cotalpha0*cotalpha0;
1625  if(arg < 0.f) arg = 0.f;
1626  cotb = std::sqrt(arg);
1627 
1628 // Copy the charge scaling factor to the private variable
1629 
1630  Ny = theStripTemp_[index_id_].head.NTy;
1631 
1632 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1633  if(Ny < 2) {
1634  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny << std::endl;
1635  }
1636 #else
1637  assert(Ny > 1);
1638 #endif
1639 
1640 // next, loop over all y-angle entries
1641 
1642  ilow = 0;
1643  yratio = 0.f;
1644 
1645  if(cotb >= theStripTemp_[index_id_].enty[Ny-1].cotbeta) {
1646 
1647  ilow = Ny-2;
1648  yratio = 1.f;
1649 
1650  } else {
1651 
1652  if(cotb >= theStripTemp_[index_id_].enty[0].cotbeta) {
1653 
1654  for (i=0; i<Ny-1; ++i) {
1655 
1656  if( theStripTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < theStripTemp_[index_id_].enty[i+1].cotbeta) {
1657 
1658  ilow = i;
1659  yratio = (cotb - theStripTemp_[index_id_].enty[i].cotbeta)/(theStripTemp_[index_id_].enty[i+1].cotbeta - theStripTemp_[index_id_].enty[i].cotbeta);
1660  break;
1661  }
1662  }
1663  }
1664  }
1665 
1666  ihigh=ilow + 1;
1667 
1668 // Interpolate Vavilov parameters
1669 
1670  mpvvav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].mpvvav + yratio*theStripTemp_[index_id_].enty[ihigh].mpvvav;
1671  sigmavav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].sigmavav + yratio*theStripTemp_[index_id_].enty[ihigh].sigmavav;
1672  kappavav_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].kappavav + yratio*theStripTemp_[index_id_].enty[ihigh].kappavav;
1673 
1674 // Copy to parameter list
1675 
1676 
1677  mpv = (double)mpvvav_;
1678  sigma = (double)sigmavav_;
1679  kappa = (double)kappavav_;
1680 
1681  return;
1682 
1683 } // vavilov_pars
1684 
1685 // ************************************************************************************************************
1690 // ************************************************************************************************************
1691 void SiStripTemplate::vavilov2_pars(double& mpv, double& sigma, double& kappa)
1692 
1693 {
1694  // Local variables
1695  int i;
1696  int ilow, ihigh, Ny;
1697  float yratio, cotb, cotalpha0, arg;
1698 
1699 // Interpolate in cotbeta only for the correct total path length (converts cotalpha, cotbeta into an effective cotbeta)
1700 
1701  cotalpha0 = theStripTemp_[index_id_].enty[0].cotalpha;
1702  arg = cotb_current_*cotb_current_ + cota_current_*cota_current_ - cotalpha0*cotalpha0;
1703  if(arg < 0.f) arg = 0.f;
1704  cotb = std::sqrt(arg);
1705 
1706  // Copy the charge scaling factor to the private variable
1707 
1708  Ny = theStripTemp_[index_id_].head.NTy;
1709 
1710 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
1711  if(Ny < 2) {
1712  throw cms::Exception("DataCorrupt") << "template ID = " << id_current_ << "has too few entries: Ny = " << Ny << std::endl;
1713  }
1714 #else
1715  assert(Ny > 1);
1716 #endif
1717 
1718  // next, loop over all y-angle entries
1719 
1720  ilow = 0;
1721  yratio = 0.f;
1722 
1723  if(cotb >= theStripTemp_[index_id_].enty[Ny-1].cotbeta) {
1724 
1725  ilow = Ny-2;
1726  yratio = 1.f;
1727 
1728  } else {
1729 
1730  if(cotb >= theStripTemp_[index_id_].enty[0].cotbeta) {
1731 
1732  for (i=0; i<Ny-1; ++i) {
1733 
1734  if( theStripTemp_[index_id_].enty[i].cotbeta <= cotb && cotb < theStripTemp_[index_id_].enty[i+1].cotbeta) {
1735 
1736  ilow = i;
1737  yratio = (cotb - theStripTemp_[index_id_].enty[i].cotbeta)/(theStripTemp_[index_id_].enty[i+1].cotbeta - theStripTemp_[index_id_].enty[i].cotbeta);
1738  break;
1739  }
1740  }
1741  }
1742  }
1743 
1744  ihigh=ilow + 1;
1745 
1746  // Interpolate Vavilov parameters
1747 
1748  mpvvav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].mpvvav2 + yratio*theStripTemp_[index_id_].enty[ihigh].mpvvav2;
1749  sigmavav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].sigmavav2 + yratio*theStripTemp_[index_id_].enty[ihigh].sigmavav2;
1750  kappavav2_ = (1.f - yratio)*theStripTemp_[index_id_].enty[ilow].kappavav2 + yratio*theStripTemp_[index_id_].enty[ihigh].kappavav2;
1751 
1752  // Copy to parameter list
1753 
1754  mpv = (double)mpvvav2_;
1755  sigma = (double)sigmavav2_;
1756  kappa = (double)kappavav2_;
1757 
1758  return;
1759 
1760 } // vavilov2_pars
1761 
1762 
1763 
SiStripTemplateEntry::clslenx
float clslenx
cluster x-length in strips at signal height sxmax/2
Definition: SiStripTemplate.h:43
SiStripTemplateEntry::xavgc2m
float xavgc2m[4]
1st pass chi2 min search: average x-bias of reconstruction binned in 4 charge bins ...
Definition: SiStripTemplate.h:61
SiStripTemplateHeader::NTy
int NTy
number of Template y entries
Definition: SiStripTemplate.h:90
i
int i
Definition: DBlmapReader.cc:9
SiStripTemplateEntry::xavggen
float xavggen[4]
generic algorithm: average x-bias of reconstruction binned in 4 charge bins
Definition: SiStripTemplate.h:67
SiStripTemplateEntry::xrmsgen
float xrmsgen[4]
generic algorithm: average x-rms of reconstruction binned in 4 charge bins
Definition: SiStripTemplate.h:68
SiStripTemplateEntry::beta
float beta
beta track angle (defined in CMS CMS IN 2004/014)
Definition: SiStripTemplate.h:34
SiStripTemplateEntry::qbfrac
float qbfrac[3]
fraction of sample in qbin = 0-2 (&gt;=3 is the complement)
Definition: SiStripTemplate.h:75
SiStripTemplateEntry::qavg_avg
float qavg_avg
average cluster charge of clusters that are less than qavg (normalize 2-D simple templates) ...
Definition: SiStripTemplate.h:77
SiStripTemplateEntry::mpvvav
float mpvvav
most probable charge in Vavilov distribution (not actually for larger kappa)
Definition: SiStripTemplate.h:44
SiStripTemplateHeader::templ_version
int templ_version
Version number of the template to ensure code compatibility.
Definition: SiStripTemplate.h:88
SiStripTemplate::xsigma2
void xsigma2(int fxstrp, int lxstrp, float sxthr, float xsum[BSXSIZE], float xsig2[BSXSIZE])
Definition: SiStripTemplate.cc:1017
SiStripTemplate::vavilov2_pars
void vavilov2_pars(double &mpv, double &sigma, double &kappa)
Definition: SiStripTemplate.cc:1691
SiStripTemplateEntry::costrk
float costrk[3]
direction cosines of tracks used to generate this entry
Definition: SiStripTemplate.h:36
SiStripTemplateHeader::temperature
float temperature
detector temperature in deg K
Definition: SiStripTemplate.h:95
SiStripTemplateEntry::xgsiggen
float xgsiggen[4]
generic algorithm: average sigma_x from Gaussian fit binned in 4 charge bins
Definition: SiStripTemplate.h:70
SiStripTemplateEntry::chi2xavgc2m
float chi2xavgc2m[4]
1st pass chi2 min search: average x chi^2 in 4 charge bins (merged clusters)
Definition: SiStripTemplate.h:65
BSXSIZE
#define BSXSIZE
Definition: SiStripTemplateDefs.h:25
SiStripTemplateStore::enty
SiStripTemplateEntry enty[31]
60 Barrel y templates spanning cluster lengths from 0px to +18px [28 entries for fstrp] ...
Definition: SiStripTemplate.h:111
SiStripTemplateHeader::ysize
float ysize
strip size (for future use in upgraded geometry)
Definition: SiStripTemplate.h:102
SiStripTemplateEntry::chi2xavgone
float chi2xavgone
average x chi^2 for 1 strip clusters
Definition: SiStripTemplate.h:59
SiStripTemplateHeader::ID
int ID
template ID number
Definition: SiStripTemplate.h:87
SiStripTemplate::xtemp3d
void xtemp3d(int j, int k, std::vector< float > &xtemplate)
Definition: SiStripTemplate.cc:1441
SiStripTemplate::vavilov_pars
void vavilov_pars(double &mpv, double &sigma, double &kappa)
Definition: SiStripTemplate.cc:1613
EcalCondDB.db
tuple db
Definition: EcalCondDB.py:151
SiStripTemplateHeader::xsize
float xsize
strip size (for future use in upgraded geometry)
Definition: SiStripTemplate.h:101
BSHX
#define BSHX
Definition: SiStripTemplateDefs.h:26
SiStripTemplateHeader::s50
float s50
1/2 of the readout threshold in ADC units
Definition: SiStripTemplate.h:98
SiStripTemplateEntry::dxone
float dxone
mean offset/correction for one strip x-clusters
Definition: SiStripTemplate.h:39
SiStripTemplateEntry::xflpar
float xflpar[4][6]
Aqfl-parameterized x-correction in 4 charge bins.
Definition: SiStripTemplate.h:56
SiStripTemplate::qbin
int qbin(int id, float cotalpha, float cotbeta, float qclus)
Definition: SiStripTemplate.cc:1469
SiStripTemplateEntry::sigmavav2
float sigmavav2
&quot;sigma&quot; scale fctor for Vavilov distribution for 2 merged clusters
Definition: SiStripTemplate.h:48
SiStripTemplateEntry::xavg
float xavg[4]
average x-bias of reconstruction binned in 4 charge bins
Definition: SiStripTemplate.h:52
SiStripTemplateEntry::runnum
int runnum
&lt; Basic template entry corresponding to a single set of track angles
Definition: SiStripTemplate.h:31
SiStripTemplateEntry::kappavav2
float kappavav2
kappa parameter for Vavilov distribution for 2 merged clusters
Definition: SiStripTemplate.h:49
ENDL
#define ENDL
Definition: SiStripTemplate.cc:38
funct::arg
A arg
Definition: Factorize.h:36
SiStripTemplateEntry::xtemp
float xtemp[9][TSXSIZE]
templates for x-reconstruction (binned over 1 central strip)
Definition: SiStripTemplate.h:51
SiStripTemplateEntry::xavgbcn
float xavgbcn[4]
barycenter: average x-bias of reconstruction binned in 4 charge bins
Definition: SiStripTemplate.h:71
SiStripTemplateHeader::zsize
float zsize
strip size (for future use in upgraded geometry)
Definition: SiStripTemplate.h:103
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
SiStripTemplateStore::entx
SiStripTemplateEntry entx[5][73]
29 Barrel x templates spanning cluster lengths from -6px (-1.125Rad) to +6px (+1.125Rad) in each of 5...
Definition: SiStripTemplate.h:112
create_public_lumi_plots.in_file
tuple in_file
Definition: create_public_lumi_plots.py:868
SiStripTemplate::xflcorr
float xflcorr(int binq, float qflx)
Definition: SiStripTemplate.cc:1132
SiStripTemplateHeader::title
char title[80]
&lt; template header structure
Definition: SiStripTemplate.h:86
SiStripTemplateEntry::chi2xavg
float chi2xavg[4]
average x chi^2 in 4 charge bins
Definition: SiStripTemplate.h:57
SiStripTemplateHeader::lorywidth
float lorywidth
estimate of y-lorentz width from single strip offset
Definition: SiStripTemplate.h:99
SiStripTemplate::pushfile
bool pushfile(int filenum)
Definition: SiStripTemplate.cc:54
SiStripTemplate::xtemp3d_int
void xtemp3d_int(int nxpix, int &nxbins)
Definition: SiStripTemplate.cc:1362
SiStripTemplateStore::head
SiStripTemplateHeader head
&lt; template storage structure
Definition: SiStripTemplate.h:109
SiStripTemplate::cxtemp
int cxtemp()
Return central pixel of x-template pixels above readout threshold.
Definition: SiStripTemplate.cc:1314
SiStripTemplateEntry::chi2xminc2m
float chi2xminc2m[4]
1st pass chi2 min search: minimum of x chi^2 in 4 charge bins (merged clusters)
Definition: SiStripTemplate.h:66
SiStripTemplateEntry::fracxone
float fracxone
fraction of sample with xsize = 1
Definition: SiStripTemplate.h:76
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
SiStripTemplateEntry::mpvvav2
float mpvvav2
most probable charge in Vavilov distribution for 2 merged clusters (not actually for larger kappa) ...
Definition: SiStripTemplate.h:47
SiStripTemplateEntry::xrmsbcn
float xrmsbcn[4]
barycenter: average x-rms of reconstruction binned in 4 charge bins
Definition: SiStripTemplate.h:72
BSXM3
#define BSXM3
Definition: SiStripTemplateDefs.h:29
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
j
int j
Definition: DBlmapReader.cc:9
SiStripTemplateEntry::sxone
float sxone
rms for one strip x-clusters
Definition: SiStripTemplate.h:40
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
LOGERROR
#define LOGERROR(x)
Definition: SiStripTemplate.cc:36
SiStripTemplateHeader::lorxwidth
float lorxwidth
estimate of x-lorentz width from single strip offset
Definition: SiStripTemplate.h:100
SiPixelTemplateDBObject::sVector
std::vector< float > sVector() const
Definition: SiPixelTemplateDBObject.h:65
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:81
SiStripTemplateHeader::qscale
float qscale
Charge scaling to match cmssw and stripav.
Definition: SiStripTemplate.h:97
SiStripTemplateEntry::xgx0c2m
float xgx0c2m[4]
1st pass chi2 min search: average x0 from Gaussian fit binned in 4 charge bins
Definition: SiStripTemplate.h:63
SiStripTemplateEntry::kappavav
float kappavav
kappa parameter for Vavilov distribution
Definition: SiStripTemplate.h:46
SiStripTemplateHeader::NTxx
int NTxx
number of Template x-entries in each slice
Definition: SiStripTemplate.h:92
SiStripTemplateEntry::chi2xminone
float chi2xminone
minimum of x chi^2 for 1 strip clusters
Definition: SiStripTemplate.h:60
SiStripTemplateEntry::qavg
float qavg
average cluster charge for this set of track angles (now includes threshold effects) ...
Definition: SiStripTemplate.h:37
SiStripTemplate::sxtemp
void sxtemp(float xhit, std::vector< float > &cluster)
Definition: SiStripTemplate.cc:1259
SiStripTemplateEntry::alpha
float alpha
alpha track angle (defined in CMS CMS IN 2004/014)
Definition: SiStripTemplate.h:32
SiStripTemplateEntry::xgx0bcn
float xgx0bcn[4]
barycenter: average x0 from Gaussian fit binned in 4 charge bins
Definition: SiStripTemplate.h:73
SiStripTemplateEntry::xgx0gen
float xgx0gen[4]
generic algorithm: average x0 from Gaussian fit binned in 4 charge bins
Definition: SiStripTemplate.h:69
SiStripTemplateEntry::xgsigbcn
float xgsigbcn[4]
barycenter: average sigma_x from Gaussian fit binned in 4 charge bins
Definition: SiStripTemplate.h:74
SiStripTemplateEntry::cotalpha
float cotalpha
cot(alpha) is proportional to cluster length in x and is basis of interpolation
Definition: SiStripTemplate.h:33
SiStripTemplateHeader::NTyx
int NTyx
number of Template y-slices of x entries
Definition: SiStripTemplate.h:91
SiStripTemplateEntry::chi2xmin
float chi2xmin[4]
minimum of x chi^2 in 4 charge bins
Definition: SiStripTemplate.h:58
SiStripTemplateEntry::sigmavav
float sigmavav
&quot;sigma&quot; scale fctor for Vavilov distribution
Definition: SiStripTemplate.h:45
edm::shift
static unsigned int const shift
Definition: LuminosityBlockID.cc:9
SiStripTemplateEntry::xgx0
float xgx0[4]
average x0 from Gaussian fit binned in 4 charge bins
Definition: SiStripTemplate.h:54
SiStripTemplateHeader::Vbias
float Vbias
detector bias potential in Volts
Definition: SiStripTemplate.h:94
SiStripTemplateEntry::xrms
float xrms[4]
average x-rms of reconstruction binned in 4 charge bins
Definition: SiStripTemplate.h:53
SiStripTemplateEntry::xpar
float xpar[2][5]
projected x-strip uncertainty parameterization
Definition: SiStripTemplate.h:50
SiStripTemplateHeader::Dtype
int Dtype
detector type (0=BPix, 1=FPix)
Definition: SiStripTemplate.h:93
LOGINFO
#define LOGINFO(x)
Definition: SiStripTemplate.cc:37
SiStripTemplate::interpolate
bool interpolate(int id, float cotalpha, float cotbeta, float locBy)
Definition: SiStripTemplate.cc:641
SiStripTemplateEntry::qmin2
float qmin2
tighter minimum cluster charge for valid hit (keeps 99.8% of simulated hits)
Definition: SiStripTemplate.h:42
SiStripTemplateEntry::sxmax
float sxmax
average strip signal for x-projection of cluster
Definition: SiStripTemplate.h:38
SiStripTemplateEntry::xrmsc2m
float xrmsc2m[4]
1st pass chi2 min search: average x-rms of reconstruction binned in 4 charge bins ...
Definition: SiStripTemplate.h:62
SiStripTemplate::xtemp
void xtemp(int fxbin, int lxbin, float xtemplate[41][BSXSIZE])
Definition: SiStripTemplate.cc:1183
TSXSIZE
#define TSXSIZE
Definition: SiStripTemplateDefs.h:22
SiStripTemplateEntry::xgsigc2m
float xgsigc2m[4]
1st pass chi2 min search: average sigma_x from Gaussian fit binned in 4 charge bins ...
Definition: SiStripTemplate.h:64
SiStripTemplateEntry::xgsig
float xgsig[4]
average sigma_x from Gaussian fit binned in 4 charge bins
Definition: SiStripTemplate.h:55
SiStripTemplateEntry::cotbeta
float cotbeta
cot(beta) is proportional to cluster length in y and is basis of interpolation
Definition: SiStripTemplate.h:35
SiStripTemplateHeader::fluence
float fluence
radiation fluence in n_eq/cm^2
Definition: SiStripTemplate.h:96
SiStripTemplateEntry::qmin
float qmin
minimum cluster charge for valid hit (keeps 99.9% of simulated hits)
Definition: SiStripTemplate.h:41
BSXM1
#define BSXM1
Definition: SiStripTemplateDefs.h:27
BSXM2
#define BSXM2
Definition: SiStripTemplateDefs.h:28
SiStripTemplateHeader::Bfield
float Bfield
Bfield in Tesla.
Definition: SiStripTemplate.h:89
Generated for CMSSW Reference Manual by doxygen 1.8.5