CMS 3D CMS Logo

EcalUncalibRecHitTimingCCAlgo.cc
Go to the documentation of this file.
2 
4  const float stopTime,
5  const float targetTimePrecision)
6  : startTime_(startTime), stopTime_(stopTime), targetTimePrecision_(targetTimePrecision) {}
7 
9  const std::vector<double>& amplitudes,
10  const EcalPedestals::Item* aped,
11  const EcalMGPAGainRatio* aGain,
12  const FullSampleVector& fullpulse,
13  EcalUncalibratedRecHit& uncalibRecHit,
14  float& errOnTime) const {
15  constexpr unsigned int nsample = EcalDataFrame::MAXSAMPLES;
16 
17  double maxamplitude = -std::numeric_limits<double>::max();
18  float pulsenorm = 0.;
19 
20  std::vector<float> pedSubSamples(nsample);
21  for (unsigned int iSample = 0; iSample < nsample; iSample++) {
22  const EcalMGPASample& sample = dataFrame.sample(iSample);
23 
24  float amplitude = 0.;
25  int gainId = sample.gainId();
26 
27  double pedestal = 0.;
28  double gainratio = 1.;
29 
30  if (gainId == 0 || gainId == 3) {
31  pedestal = aped->mean_x1;
32  gainratio = aGain->gain6Over1() * aGain->gain12Over6();
33  } else if (gainId == 1) {
34  pedestal = aped->mean_x12;
35  gainratio = 1.;
36  } else if (gainId == 2) {
37  pedestal = aped->mean_x6;
38  gainratio = aGain->gain12Over6();
39  }
40 
41  amplitude = (static_cast<float>(sample.adc()) - pedestal) * gainratio;
42 
43  if (gainId == 0) {
44  //saturation
45  amplitude = (4095. - pedestal) * gainratio;
46  }
47 
48  pedSubSamples[iSample] = amplitude;
49 
50  if (amplitude > maxamplitude) {
51  maxamplitude = amplitude;
52  }
53  pulsenorm += fullpulse(iSample);
54  }
55 
56  int ipulse = -1;
57  for (auto const& amplit : amplitudes) {
58  ipulse++;
59  int bxp3 = ipulse - 2;
60  int firstsamplet = std::max(0, bxp3);
61  int offset = 7 - bxp3;
62 
63  for (unsigned int isample = firstsamplet; isample < nsample; ++isample) {
64  auto const pulse = fullpulse(isample + offset);
65  pedSubSamples[isample] = std::max(0., pedSubSamples[isample] - amplit * pulse / pulsenorm);
66  }
67  }
68 
69  // Start of time computation
70  float tStart = startTime_ + GLOBAL_TIME_SHIFT;
71  float tStop = stopTime_ + GLOBAL_TIME_SHIFT;
72  float tM = (tStart + tStop) / 2;
73 
74  float distStart, distStop;
75  int counter = 0;
76 
77  do {
78  ++counter;
79  distStart = computeCC(pedSubSamples, fullpulse, tStart);
80  distStop = computeCC(pedSubSamples, fullpulse, tStop);
81 
82  if (distStart > distStop) {
83  tStop = tM;
84  } else {
85  tStart = tM;
86  }
87  tM = (tStart + tStop) / 2;
88 
89  } while (tStop - tStart > targetTimePrecision_ && counter < MAX_NUM_OF_ITERATIONS);
90 
91  tM -= GLOBAL_TIME_SHIFT;
92  errOnTime = targetTimePrecision_;
93 
94  if (counter < MIN_NUM_OF_ITERATIONS || counter > MAX_NUM_OF_ITERATIONS - 1) {
96  //Negative error means that there was a problem with the CC
98  }
99 
100  return -tM / ecalPh1::Samp_Period;
101 }
102 
104  const float time) const {
105  // t is in ns
107  if (time < 0)
108  shift -= 1;
109  float tt = time / ecalPh1::Samp_Period - shift;
110 
111  FullSampleVector interpPulse;
112  // 2nd poly with avg
113  unsigned int numberOfSamples = fullpulse.size();
114  auto facM1orP2 = 0.25 * tt * (tt - 1);
115  auto fac = (0.25 * (tt - 2) - 0.5 * (tt + 1)) * (tt - 1);
116  auto facP1 = (0.25 * (tt + 1) - 0.5 * (tt - 2)) * tt;
117  for (unsigned int i = 1; i < numberOfSamples - 2; ++i) {
118  float a =
119  facM1orP2 * fullpulse[i - 1] + fac * fullpulse[i] + facP1 * fullpulse[i + 1] + facM1orP2 * fullpulse[i + 2];
120  if (a > 0)
121  interpPulse[i] = a;
122  else
123  interpPulse[i] = 0;
124  }
125  interpPulse[0] = facM1orP2 * fullpulse[0] + facP1 * fullpulse[1] + facM1orP2 * fullpulse[2];
126  interpPulse[numberOfSamples - 2] = facM1orP2 * fullpulse[numberOfSamples - 3] + fac * fullpulse[numberOfSamples - 2] +
127  facP1 * fullpulse[numberOfSamples - 1];
128  interpPulse[numberOfSamples - 1] = 2 * facM1orP2 * fullpulse[numberOfSamples - 2] -
129  4 * facM1orP2 * fullpulse[numberOfSamples - 1] +
130  facP1 * fullpulse[numberOfSamples - 1];
131 
132  FullSampleVector interpPulseShifted;
133  for (int i = 0; i < interpPulseShifted.size(); ++i) {
134  if (i + shift >= 0 && i + shift < interpPulse.size())
135  interpPulseShifted[i] = interpPulse[i + shift];
136  else
137  interpPulseShifted[i] = 0;
138  }
139  return interpPulseShifted;
140 }
141 
142 float EcalUncalibRecHitTimingCCAlgo::computeCC(const std::vector<float>& samples,
143  const FullSampleVector& signalTemplate,
144  const float time) const {
145  constexpr int exclude = 1;
146  float powerSamples = 0.;
147  float powerTemplate = 0.;
148  float cc = 0.;
149  auto interpolated = interpolatePulse(signalTemplate, time);
150  for (int i = exclude; i < int(samples.size() - exclude); ++i) {
151  powerSamples += std::pow(samples[i], 2);
152  powerTemplate += std::pow(interpolated[i], 2);
153  cc += interpolated[i] * samples[i];
154  }
155 
156  float denominator = std::sqrt(powerTemplate * powerSamples);
157  return cc / denominator;
158 }
Eigen::Matrix< double, FullSampleVectorSize, 1 > FullSampleVector
string startTime
Definition: TTTypes.h:54
T sqrt(T t)
Definition: SSEVec.h:19
EcalMGPASample sample(int i) const
Definition: EcalDataFrame.h:29
static constexpr double Samp_Period
Definition: EcalConstants.h:41
float computeCC(const std::vector< float > &samples, const FullSampleVector &sigmalTemplate, const float t) const
float gain12Over6() const
constexpr int gainId(sample_type sample)
get the gainId (2 bits)
double pulse(double x, double y, double z, double t)
float gain6Over1() const
double a
Definition: hdecay.h:119
static std::atomic< unsigned int > counter
static unsigned int const shift
double computeTimeCC(const EcalDataFrame &dataFrame, const std::vector< double > &amplitudes, const EcalPedestals::Item *aped, const EcalMGPAGainRatio *aGain, const FullSampleVector &fullpulse, EcalUncalibratedRecHit &uncalibRecHit, float &errOnTime) const
static constexpr int MAXSAMPLES
Definition: EcalDataFrame.h:48
FullSampleVector interpolatePulse(const FullSampleVector &fullpulse, const float t=0) const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
EcalUncalibRecHitTimingCCAlgo(const float startTime, const float stopTime, const float targetTimePrecision)