CMS 3D CMS Logo

MSLayersKeeperX0AtEta.cc
Go to the documentation of this file.
4 #include <algorithm>
6 
7 using namespace std;
8 
9 namespace {
10  template <class T>
11  inline T sqr(T t) {
12  return t * t;
13  }
14  // float unsafe_asinhf(float x) { return std::abs(x)>10.f ? std::copysign(3.f,x) : unsafe_logf<3>(x+std::sqrt(1.f+x*x)); }
15  inline float unsafe_asinhf(float x) { return unsafe_logf<3>(x + std::sqrt(1.f + x * x)); }
16 } // namespace
17 
18 //------------------------------------------------------------------------------
19 const MSLayersAtAngle& MSLayersKeeperX0AtEta::layers(float cotTheta) const {
20  float eta = unsafe_asinhf(cotTheta);
21  return theLayersData[idxBin(eta)];
22 }
23 
24 //------------------------------------------------------------------------------
25 float MSLayersKeeperX0AtEta::eta(int idxBin) const { return (idxBin + 0.5 - theHalfNBins) * theDeltaEta; }
26 
27 //------------------------------------------------------------------------------
29  float ieta = eta / theDeltaEta;
30  if (std::abs(ieta) >= theHalfNBins - 1.e-3)
31  return (eta > 0) ? max(2 * theHalfNBins - 1, 0) : 0;
32  else
33  return int(ieta + theHalfNBins);
34 }
35 
36 //------------------------------------------------------------------------------
38  if (isInitialised)
39  return;
40  isInitialised = true;
41  const float BIG = 99999.;
42 
43  // set size from data file
44  MultipleScatteringX0Data msX0data;
45  theHalfNBins = msX0data.nBinsEta();
46  float etaMax = msX0data.maxEta();
47  theDeltaEta = (theHalfNBins != 0) ? etaMax / theHalfNBins : BIG;
48 
49  theLayersData = vector<MSLayersAtAngle>(max(2 * theHalfNBins, 1));
51  for (int idxbin = 0; idxbin < 2 * theHalfNBins; idxbin++) {
52  float etaValue = eta(idxbin);
53  float cotTheta = sinh(etaValue);
54 
55  vector<MSLayer> layers = layout.detLayers(etaValue, 0, iSetup);
56  vector<MSLayer> tmplay = layout.otherLayers(etaValue, iSetup);
57  layers.insert(layers.end(), tmplay.begin(), tmplay.end());
58  sort(layers.begin(), layers.end());
59  setX0(layers, etaValue, msX0data);
60  theLayersData[idxbin] = MSLayersAtAngle(layers);
61  PixelRecoPointRZ zero(0., 0.);
62  PixelRecoLineRZ line(zero, cotTheta);
63  vector<MSLayer>::iterator it;
64  for (it = layers.begin(); it != layers.end(); it++) {
65  float x0 = getDataX0(*it).x0;
66  float sumX0D = theLayersData[idxbin].sumX0D(zero, it->crossing(line).first);
67  setDataX0(*it, DataX0(x0, sumX0D, cotTheta));
68  theLayersData[idxbin].update(*it);
69  }
70  }
71 
72  // add layers not seen from nominal vertex but crossed if
73  // vertex seperated from nominal by less than 3 sigma
74  for (int idxbin = 0; idxbin < 2 * theHalfNBins; idxbin++) {
75  float etaValue = eta(idxbin);
76  for (int isign = -1; isign <= 1; isign += 2) {
77  float z = isign * 15.9; //3 sigma from zero
78  const MSLayersAtAngle& layersAtAngle = theLayersData[idxbin];
79  vector<MSLayer> candidates = layout.detLayers(etaValue, z, iSetup);
80  vector<MSLayer>::iterator it;
81  for (it = candidates.begin(); it != candidates.end(); it++) {
82  if (layersAtAngle.findLayer(*it))
83  continue;
84  const MSLayer* found = nullptr;
85  int bin = idxbin;
86  while (!found) {
87  bin--;
88  if (bin < 0)
89  break;
90  found = theLayersData[bin].findLayer(*it);
91  }
92  bin = idxbin;
93  while (!found) {
94  bin++;
95  if (bin > 2 * theHalfNBins - 1)
96  break;
97  found = theLayersData[bin].findLayer(*it);
98  }
99  if (found)
100  theLayersData[idxbin].update(*found);
101  }
102  }
103  }
104 
105  // cout << "LAYERS, size=: "<<theLayersData.size()<< endl;
106  /*
107  for (int idxbin = 0; idxbin <= theHalfNBins; idxbin+=25) {
108  float etaValue = eta(idxbin);
109  const MSLayersAtAngle & layers= theLayersData[idxbin];
110  cout << "ETA: "<< etaValue <<" (bin:"<<idxbin<<") #layers:"
111  <<layers.size()<<endl;
112  layers.print();
113  }
114  for (int idxbin = 2*theHalfNBins-1; idxbin > theHalfNBins; idxbin-=25) {
115  float etaValue = eta(idxbin);
116  const MSLayersAtAngle & layers= theLayersData[idxbin];
117  cout << "ETA: "<< etaValue <<" (bin:"<<idxbin<<") #layers:"
118  <<layers.size()<<endl;
119  layers.print();
120  }
121 */
122 }
123 
124 //------------------------------------------------------------------------------
125 void MSLayersKeeperX0AtEta::setX0(vector<MSLayer>& layers, float eta, const SumX0AtEtaDataProvider& sumX0) {
126  const float BIG = 99999.;
127  float cotTheta = sinh(eta);
128  float sinTheta = 1 / sqrt(1 + sqr(cotTheta));
129  float cosTheta = cotTheta * sinTheta;
130 
131  float sumX0atAngleLast = 0.;
132  vector<MSLayer>::iterator il;
133  for (il = layers.begin(); il != layers.end(); il++) {
134  PixelRecoLineRZ line(PixelRecoPointRZ(0., 0.), cotTheta);
135  float rL = (*il).crossing(line).first.r();
136  float rN = (il + 1 != layers.end()) ? (il + 1)->crossing(line).first.r() : BIG;
137  float rBound = (rL + rN) / 2.;
138  float sumX0atAngle = sumX0.sumX0atEta(eta, rBound);
139 
140  float dX0 = (il->face() == GeomDetEnumerators::barrel) ? (sumX0atAngle - sumX0atAngleLast) * sinTheta
141  : (sumX0atAngle - sumX0atAngleLast) * fabs(cosTheta);
142 
143  setDataX0(*il, DataX0(dX0, 0, cotTheta));
144  sumX0atAngleLast = sumX0atAngle;
145  }
146 }
SumX0AtEtaDataProvider
Definition: MultipleScatteringX0Data.h:13
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
MSLayersKeeperX0AtEta::eta
float eta(int idxBin) const
Definition: MSLayersKeeperX0AtEta.cc:25
sqr
int sqr(const T &t)
Definition: pfalgo_common_ref.h:9
MultipleScatteringGeometry
Definition: MultipleScatteringGeometry.h:11
makeLayoutFileForGui.layout
layout
Definition: makeLayoutFileForGui.py:67
SiPixelPI::zero
Definition: SiPixelPayloadInspectorHelper.h:39
MultipleScatteringX0Data::maxEta
float maxEta() const
Definition: MultipleScatteringX0Data.cc:50
newFWLiteAna.found
found
Definition: newFWLiteAna.py:118
GeomDetEnumerators::barrel
Definition: GeomDetEnumerators.h:9
MSLayersKeeperX0AtEta::layers
const MSLayersAtAngle & layers(float cotTheta) const override
Definition: MSLayersKeeperX0AtEta.cc:19
MSLayersAtAngle::findLayer
const MSLayer * findLayer(const MSLayer &layer) const
Definition: MSLayersAtAngle.cc:34
MultipleScatteringX0Data::nBinsEta
int nBinsEta() const
Definition: MultipleScatteringX0Data.cc:36
MultipleScatteringGeometry.h
PVValHelper::eta
Definition: PVValidationHelpers.h:70
MSLayersKeeperX0AtEta.h
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
MultipleScatteringX0Data
Definition: MultipleScatteringX0Data.h:19
MSLayersAtAngle
Definition: MSLayersAtAngle.h:16
MSLayer::DataX0
Definition: MSLayer.h:17
MSLayer
Definition: MSLayer.h:13
approx_log.h
MSLayersKeeperX0AtEta::idxBin
int idxBin(float eta) const
Definition: MSLayersKeeperX0AtEta.cc:28
MSLayersKeeperX0AtEta::init
void init(const edm::EventSetup &iSetup) override
Definition: MSLayersKeeperX0AtEta.cc:37
LEDCalibrationChannels.ieta
ieta
Definition: LEDCalibrationChannels.py:63
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
createfilelist.int
int
Definition: createfilelist.py:10
PixelRecoLineRZ
Definition: PixelRecoLineRZ.h:12
edm::EventSetup
Definition: EventSetup.h:58
PixelRecoPointRZ
Definition: PixelRecoPointRZ.h:6
newFWLiteAna.bin
bin
Definition: newFWLiteAna.py:161
std
Definition: JetResolutionObject.h:76
T
long double T
Definition: Basic3DVectorLD.h:48
HLT_FULL_cff.candidates
candidates
Definition: HLT_FULL_cff.py:54992
SumX0AtEtaDataProvider::sumX0atEta
virtual float sumX0atEta(float eta, float r) const =0
MultipleScatteringX0Data.h
ALCARECOTkAlBeamHalo_cff.etaMax
etaMax
Definition: ALCARECOTkAlBeamHalo_cff.py:33
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
MSLayersKeeperX0AtEta::setX0
static void setX0(std::vector< MSLayer > &, float eta, const SumX0AtEtaDataProvider &)
Definition: MSLayersKeeperX0AtEta.cc:125
mps_splice.line
line
Definition: mps_splice.py:76
submitPVValidationJobs.t
string t
Definition: submitPVValidationJobs.py:644
hgcalTopologyTester_cfi.layers
layers
Definition: hgcalTopologyTester_cfi.py:8
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37