CMS 3D CMS Logo
CalibrationAlgorithm.cc
Go to the documentation of this file.
1 #include "DQM/SiStripCommissioningAnalysis/interface/CalibrationAlgorithm.h"
2 #include "CondFormats/SiStripObjects/interface/CalibrationAnalysis.h"
3 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
4 #include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h"
5 #include "FWCore/MessageLogger/interface/MessageLogger.h"
6 #include "DQM/SiStripCommissioningAnalysis/interface/SiStripPulseShape.h"
7 #include "TProfile.h"
8 #include "TF1.h"
9 #include "TH1.h"
10 #include "TVirtualFitter.h"
11 #include "TFitResultPtr.h"
12 #include "TFitResult.h"
13 #include <iostream>
14 #include <sstream>
15 #include <iomanip>
16 #include <cmath>
17 #include "Math/MinimizerOptions.h"
18 
19 using namespace sistrip;
20 
21 // ----------------------------------------------------------------------------
22 //
23 CalibrationAlgorithm::CalibrationAlgorithm(const edm::ParameterSet& pset, CalibrationAnalysis* const anal)
24  : CommissioningAlgorithm(anal), cal_(nullptr) {}
25 
26 // ----------------------------------------------------------------------------
27 //
28 void CalibrationAlgorithm::extract(const std::vector<TH1*>& histos) {
29  // extract analysis object which should be already created
30  if (!anal()) {
31  edm::LogWarning(mlCommissioning_) << "[CalibrationAlgorithm::" << __func__ << "]"
32  << " NULL pointer to base Analysis object!";
33  return;
34  }
35 
36  CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
37  cal_ = dynamic_cast<CalibrationAnalysis*>(tmp);
38 
39  if (!cal_) {
40  edm::LogWarning(mlCommissioning_) << "[CalibrationAlgorithm::" << __func__ << "]"
41  << " NULL pointer to derived Analysis object!";
42  return;
43  }
44 
45  // Extract FED key from histo title
46  if (!histos.empty()) {
47  cal_->fedKey(extractFedKey(histos.front()));
48  }
49 
50  // Extract histograms
51  std::vector<TH1*>::const_iterator ihis = histos.begin();
52  unsigned int cnt = 0;
53  for (; ihis != histos.end(); ihis++, cnt++) {
54  // Check for NULL pointer
55  if (!(*ihis)) {
56  continue;
57  }
58 
59  // Check name
60  SiStripHistoTitle title((*ihis)->GetName());
61  if (title.runType() != sistrip::CALIBRATION && title.runType() != sistrip::CALIBRATION_DECO) {
62  cal_->addErrorCode(sistrip::unexpectedTask_);
63  continue;
64  }
65 
67  std::vector<std::string> tokens;
68  std::string token;
69  std::istringstream tokenStream(title.extraInfo());
70  while (std::getline(tokenStream, token, '_')) {
71  tokens.push_back(token);
72  }
73 
75  Histo histo_temp;
76  histo_temp.first = *ihis;
77  histo_temp.second = (*ihis)->GetTitle();
78  histo_temp.first->Sumw2();
79  histo_.push_back(histo_temp);
80  apvId_.push_back(title.channel() % 2);
81  stripId_.push_back(std::stoi(tokens.at(1)) * 16 + std::stoi(tokens.at(3)));
82  calChan_.push_back(std::stoi(tokens.at(1)));
83  }
84 }
85 
86 // ----------------------------------------------------------------------------
87 //
88 void CalibrationAlgorithm::analyse() {
89  ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
90  ROOT::Math::MinimizerOptions::SetDefaultStrategy(0);
91 
92  if (!cal_) {
93  edm::LogWarning(mlCommissioning_) << "[CalibrationAlgorithm::" << __func__ << "]"
94  << " NULL pointer to derived Analysis object!";
95  return;
96  }
97 
98  float Amean[2] = {-1., -1.};
99  float Amin[2] = {-1., -1.};
100  float Amax[2] = {-1., -1.};
101  float Aspread[2] = {-1., -1.};
102  float Tmean[2] = {-1., -1.};
103  float Tmin[2] = {-1., -1.};
104  float Tmax[2] = {-1., -1.};
105  float Tspread[2] = {-1., -1.};
106  float Rmean[2] = {-1., -1.};
107  float Rmin[2] = {-1., -1.};
108  float Rmax[2] = {-1., -1.};
109  float Rspread[2] = {-1., -1.};
110  float Cmean[2] = {-1., -1.};
111  float Cmin[2] = {-1., -1.};
112  float Cmax[2] = {-1., -1.};
113  float Cspread[2] = {-1., -1.};
114  float Smean[2] = {-1., -1.};
115  float Smin[2] = {-1., -1.};
116  float Smax[2] = {-1., -1.};
117  float Sspread[2] = {-1., -1.};
118  float Kmean[2] = {-1., -1.};
119  float Kmin[2] = {-1., -1.};
120  float Kmax[2] = {-1., -1.};
121  float Kspread[2] = {-1., -1.};
122  // turnOn
123  float Omean[2] = {-1., -1.};
124  float Omin[2] = {-1., -1.};
125  float Omax[2] = {-1., -1.};
126  float Ospread[2] = {-1., -1.};
127  // maximum
128  float Mmean[2] = {-1., -1.};
129  float Mmin[2] = {-1., -1.};
130  float Mmax[2] = {-1., -1.};
131  float Mspread[2] = {-1., -1.};
132  // undershoot
133  float Umean[2] = {-1., -1.};
134  float Umin[2] = {-1., -1.};
135  float Umax[2] = {-1., -1.};
136  float Uspread[2] = {-1., -1.};
137  // baseline
138  float Bmean[2] = {-1., -1.};
139  float Bmin[2] = {-1., -1.};
140  float Bmax[2] = {-1., -1.};
141  float Bspread[2] = {-1., -1.};
142 
144  TFitResultPtr fit_result;
145  TF1* fit_function = nullptr;
146  if (cal_->deconv_) {
147  fit_function = new TF1("fit_function_deco", fdeconv, 0, 400, 7);
148  fit_function->SetParameters(4, 25, 25, 50, 250, 25, 0.75);
149  } else {
150  fit_function = new TF1("fit_function_peak", fpeak, 0, 400, 6);
151  fit_function->SetParameters(4, 50, 50, 70, 250, 20);
152  }
153 
155  std::vector<unsigned int> nStrips(2, 0.);
156 
157  for (size_t ihist = 0; ihist < histo_.size(); ihist++) {
158  if (!histo_[ihist].first) {
159  edm::LogWarning(mlCommissioning_) << " NULL pointer to histogram for: " << histo_[ihist].second << " !";
160  return;
161  }
162 
163  cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] = 0;
164  cal_->baseline_[apvId_[ihist]][stripId_[ihist]] = 0;
165  cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] = 0;
166  cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] = 0;
167  cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] = 0;
168  cal_->undershoot_[apvId_[ihist]][stripId_[ihist]] = 0;
169  cal_->tail_[apvId_[ihist]][stripId_[ihist]] = 0;
170  cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] = 0;
171  cal_->smearing_[apvId_[ihist]][stripId_[ihist]] = 0;
172  cal_->chi2_[apvId_[ihist]][stripId_[ihist]] = 0;
173  cal_->isvalid_[apvId_[ihist]][stripId_[ihist]] = true;
174 
175  if (histo_[ihist].first->Integral() == 0) {
176  cal_->isvalid_[apvId_[ihist]][stripId_[ihist]] = false;
177  continue;
178  }
179 
180  // rescale the plot and set reasonable errors
181  correctDistribution(histo_[ihist].first);
182 
183  // from NOTE2009_021 : The charge injection provided by the calibration circuit is known with a precision of 5%
184  float error = histo_[ihist].first->GetMaximum() * 0.05;
185  for (int i = 1; i <= histo_[ihist].first->GetNbinsX(); ++i)
186  histo_[ihist].first->SetBinError(i, error);
187 
188  // set intial par
189  if (cal_->deconv_)
190  fit_function->SetParameters(10, 15, 30, 10, 350, 50, 0.75);
191  else
192  fit_function->SetParameters(6, 40, 40, 70, 350, 20);
193 
194  fit_result = histo_[ihist].first->Fit(fit_function, "QS");
195 
196  // fit-result should exist and have a resonably good status
197  if (not fit_result.Get())
198  continue;
199 
200  float maximum_ampl = fit_function->GetMaximum();
201  float peak_time = fit_function->GetMaximumX();
202  float baseline = baseLine(fit_function);
203  float turn_on_time = turnOn(fit_function, baseline);
204  float rise_time = peak_time - turn_on_time;
205 
206  // start filling info
207  cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] = maximum_ampl - baseline;
208  cal_->baseline_[apvId_[ihist]][stripId_[ihist]] = baseline;
209  cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] = rise_time;
210  cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] = turn_on_time;
211  cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] = peak_time;
212 
213  if (cal_->deconv_ and fit_function->GetMinimumX() > peak_time) // make sure the minimum is after the peak-time
214  cal_->undershoot_[apvId_[ihist]][stripId_[ihist]] =
215  100 * (fit_function->GetMinimum() - baseline) / (maximum_ampl - baseline);
216  else
217  cal_->undershoot_[apvId_[ihist]][stripId_[ihist]] = 0;
218 
219  // Bin related to peak + 125 ns
220  int lastBin = histo_[ihist].first->FindBin(peak_time + 125);
221  if (lastBin > histo_[ihist].first->GetNbinsX() - 4)
222  lastBin = histo_[ihist].first->GetNbinsX() - 4;
223 
224  // tail is the amplitude at 5 bx from the maximum
225  cal_->tail_[apvId_[ihist]][stripId_[ihist]] =
226  100 * (histo_[ihist].first->GetBinContent(lastBin) - baseline) / (maximum_ampl - baseline);
227 
228  // reaches 1/e of the peak amplitude
229  cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] = decayTime(fit_function) - peak_time;
230  cal_->smearing_[apvId_[ihist]][stripId_[ihist]] = 0;
231  cal_->chi2_[apvId_[ihist]][stripId_[ihist]] =
232  fit_function->GetChisquare() / (histo_[ihist].first->GetNbinsX() - fit_function->GetNpar());
233 
234  // calibration channel
235  cal_->calChan_ = calChan_[ihist];
236 
237  // apply quality requirements
238  bool isvalid = true;
239  if (not cal_->deconv_) { // peak-mode
240 
241  if (cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minAmplitudeThreshold_)
242  isvalid = false;
243  if (cal_->baseline_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minBaselineThreshold_)
244  isvalid = false;
245  else if (cal_->baseline_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxBaselineThreshold_)
246  isvalid = false;
247  if (cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minDecayTimeThreshold_)
248  isvalid = false;
249  else if (cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxDecayTimeThreshold_)
250  isvalid = false;
251  if (cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minPeakTimeThreshold_)
252  isvalid = false;
253  else if (cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxPeakTimeThreshold_)
254  isvalid = false;
255  if (cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minRiseTimeThreshold_)
256  isvalid = false;
257  else if (cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxRiseTimeThreshold_)
258  isvalid = false;
259  if (cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minTurnOnThreshold_)
260  isvalid = false;
261  else if (cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxTurnOnThreshold_)
262  isvalid = false;
263  if (cal_->chi2_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxChi2Threshold_)
264  isvalid = false;
265 
266  } else {
267  if (fit_function->GetMinimumX() < peak_time)
268  isvalid = false;
269  if (cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minAmplitudeThreshold_)
270  isvalid = false;
271  if (cal_->baseline_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minBaselineThreshold_)
272  isvalid = false;
273  if (cal_->baseline_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minBaselineThreshold_)
274  isvalid = false;
275  else if (cal_->baseline_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxBaselineThreshold_)
276  isvalid = false;
277  if (cal_->chi2_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxChi2Threshold_)
278  isvalid = false;
279  if (cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minTurnOnThresholdDeco_)
280  isvalid = false;
281  else if (cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxTurnOnThresholdDeco_)
282  isvalid = false;
283  if (cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minDecayTimeThresholdDeco_)
284  isvalid = false;
285  else if (cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxDecayTimeThresholdDeco_)
286  isvalid = false;
287  if (cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minPeakTimeThresholdDeco_)
288  isvalid = false;
289  else if (cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxPeakTimeThresholdDeco_)
290  isvalid = false;
291  if (cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minRiseTimeThresholdDeco_)
292  isvalid = false;
293  else if (cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxRiseTimeThresholdDeco_)
294  isvalid = false;
295  }
296 
297  if (not isvalid) { // not valid set default to zero for all quantities
298 
299  cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] = 0;
300  cal_->baseline_[apvId_[ihist]][stripId_[ihist]] = 0;
301  cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] = 0;
302  cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] = 0;
303  cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] = 0;
304  cal_->undershoot_[apvId_[ihist]][stripId_[ihist]] = 0;
305  cal_->tail_[apvId_[ihist]][stripId_[ihist]] = 0;
306  cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] = 0;
307  cal_->smearing_[apvId_[ihist]][stripId_[ihist]] = 0;
308  cal_->chi2_[apvId_[ihist]][stripId_[ihist]] = 0;
309  cal_->isvalid_[apvId_[ihist]][stripId_[ihist]] = false;
310  continue;
311  }
312 
313  // in case is valid
314  nStrips[apvId_[ihist]]++;
315 
316  //compute mean, max, min, spread only for valid strips
317  Amean[apvId_[ihist]] += cal_->amplitude_[apvId_[ihist]][stripId_[ihist]];
318  Amin[apvId_[ihist]] = Amin[apvId_[ihist]] < cal_->amplitude_[apvId_[ihist]][stripId_[ihist]]
319  ? Amin[apvId_[ihist]]
320  : cal_->amplitude_[apvId_[ihist]][stripId_[ihist]];
321  Amax[apvId_[ihist]] = Amax[apvId_[ihist]] > cal_->amplitude_[apvId_[ihist]][stripId_[ihist]]
322  ? Amax[apvId_[ihist]]
323  : cal_->amplitude_[apvId_[ihist]][stripId_[ihist]];
324  Aspread[apvId_[ihist]] +=
325  cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] * cal_->amplitude_[apvId_[ihist]][stripId_[ihist]];
326 
327  Tmean[apvId_[ihist]] += cal_->tail_[apvId_[ihist]][stripId_[ihist]];
328  Tmin[apvId_[ihist]] = Tmin[apvId_[ihist]] < cal_->tail_[apvId_[ihist]][stripId_[ihist]]
329  ? Tmin[apvId_[ihist]]
330  : cal_->tail_[apvId_[ihist]][stripId_[ihist]];
331  Tmax[apvId_[ihist]] = Tmax[apvId_[ihist]] > cal_->tail_[apvId_[ihist]][stripId_[ihist]]
332  ? Tmax[apvId_[ihist]]
333  : cal_->tail_[apvId_[ihist]][stripId_[ihist]];
334  Tspread[apvId_[ihist]] += cal_->tail_[apvId_[ihist]][stripId_[ihist]] * cal_->tail_[apvId_[ihist]][stripId_[ihist]];
335 
336  Rmean[apvId_[ihist]] += cal_->riseTime_[apvId_[ihist]][stripId_[ihist]];
337  Rmin[apvId_[ihist]] = Rmin[apvId_[ihist]] < cal_->riseTime_[apvId_[ihist]][stripId_[ihist]]
338  ? Rmin[apvId_[ihist]]
339  : cal_->riseTime_[apvId_[ihist]][stripId_[ihist]];
340  Rmax[apvId_[ihist]] = Rmax[apvId_[ihist]] > cal_->riseTime_[apvId_[ihist]][stripId_[ihist]]
341  ? Rmax[apvId_[ihist]]
342  : cal_->riseTime_[apvId_[ihist]][stripId_[ihist]];
343  Rspread[apvId_[ihist]] +=
344  cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] * cal_->riseTime_[apvId_[ihist]][stripId_[ihist]];
345 
346  Cmean[apvId_[ihist]] += cal_->decayTime_[apvId_[ihist]][stripId_[ihist]];
347  Cmin[apvId_[ihist]] = Cmin[apvId_[ihist]] < cal_->decayTime_[apvId_[ihist]][stripId_[ihist]]
348  ? Cmin[apvId_[ihist]]
349  : cal_->decayTime_[apvId_[ihist]][stripId_[ihist]];
350  Cmax[apvId_[ihist]] = Cmax[apvId_[ihist]] > cal_->decayTime_[apvId_[ihist]][stripId_[ihist]]
351  ? Cmax[apvId_[ihist]]
352  : cal_->decayTime_[apvId_[ihist]][stripId_[ihist]];
353  Cspread[apvId_[ihist]] +=
354  cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] * cal_->decayTime_[apvId_[ihist]][stripId_[ihist]];
355 
356  Smean[apvId_[ihist]] += cal_->smearing_[apvId_[ihist]][stripId_[ihist]];
357  Smin[apvId_[ihist]] = Smin[apvId_[ihist]] < cal_->smearing_[apvId_[ihist]][stripId_[ihist]]
358  ? Smin[apvId_[ihist]]
359  : cal_->smearing_[apvId_[ihist]][stripId_[ihist]];
360  Smax[apvId_[ihist]] = Smax[apvId_[ihist]] > cal_->smearing_[apvId_[ihist]][stripId_[ihist]]
361  ? Smax[apvId_[ihist]]
362  : cal_->smearing_[apvId_[ihist]][stripId_[ihist]];
363  Sspread[apvId_[ihist]] +=
364  cal_->smearing_[apvId_[ihist]][stripId_[ihist]] * cal_->smearing_[apvId_[ihist]][stripId_[ihist]];
365 
366  Kmean[apvId_[ihist]] += cal_->chi2_[apvId_[ihist]][stripId_[ihist]];
367  Kmin[apvId_[ihist]] = Kmin[apvId_[ihist]] < cal_->chi2_[apvId_[ihist]][stripId_[ihist]]
368  ? Kmin[apvId_[ihist]]
369  : cal_->chi2_[apvId_[ihist]][stripId_[ihist]];
370  Kmax[apvId_[ihist]] = Kmax[apvId_[ihist]] > cal_->chi2_[apvId_[ihist]][stripId_[ihist]]
371  ? Kmax[apvId_[ihist]]
372  : cal_->chi2_[apvId_[ihist]][stripId_[ihist]];
373  Kspread[apvId_[ihist]] += cal_->chi2_[apvId_[ihist]][stripId_[ihist]] * cal_->chi2_[apvId_[ihist]][stripId_[ihist]];
374 
375  Omean[apvId_[ihist]] += cal_->turnOn_[apvId_[ihist]][stripId_[ihist]];
376  Omin[apvId_[ihist]] = Omin[apvId_[ihist]] < cal_->turnOn_[apvId_[ihist]][stripId_[ihist]]
377  ? Omin[apvId_[ihist]]
378  : cal_->turnOn_[apvId_[ihist]][stripId_[ihist]];
379  Omax[apvId_[ihist]] = Omax[apvId_[ihist]] > cal_->turnOn_[apvId_[ihist]][stripId_[ihist]]
380  ? Omax[apvId_[ihist]]
381  : cal_->turnOn_[apvId_[ihist]][stripId_[ihist]];
382  Ospread[apvId_[ihist]] +=
383  cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] * cal_->turnOn_[apvId_[ihist]][stripId_[ihist]];
384 
385  Mmean[apvId_[ihist]] += cal_->peakTime_[apvId_[ihist]][stripId_[ihist]];
386  Mmin[apvId_[ihist]] = Mmin[apvId_[ihist]] < cal_->peakTime_[apvId_[ihist]][stripId_[ihist]]
387  ? Mmin[apvId_[ihist]]
388  : cal_->peakTime_[apvId_[ihist]][stripId_[ihist]];
389  Mmax[apvId_[ihist]] = Mmax[apvId_[ihist]] > cal_->peakTime_[apvId_[ihist]][stripId_[ihist]]
390  ? Mmax[apvId_[ihist]]
391  : cal_->peakTime_[apvId_[ihist]][stripId_[ihist]];
392  Mspread[apvId_[ihist]] +=
393  cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] * cal_->peakTime_[apvId_[ihist]][stripId_[ihist]];
394 
395  Umean[apvId_[ihist]] += cal_->undershoot_[apvId_[ihist]][stripId_[ihist]];
396  Umin[apvId_[ihist]] = Umin[apvId_[ihist]] < cal_->undershoot_[apvId_[ihist]][stripId_[ihist]]
397  ? Umin[apvId_[ihist]]
398  : cal_->undershoot_[apvId_[ihist]][stripId_[ihist]];
399  Umax[apvId_[ihist]] = Umax[apvId_[ihist]] > cal_->undershoot_[apvId_[ihist]][stripId_[ihist]]
400  ? Umax[apvId_[ihist]]
401  : cal_->undershoot_[apvId_[ihist]][stripId_[ihist]];
402  Uspread[apvId_[ihist]] +=
403  cal_->undershoot_[apvId_[ihist]][stripId_[ihist]] * cal_->undershoot_[apvId_[ihist]][stripId_[ihist]];
404 
405  Bmean[apvId_[ihist]] += cal_->baseline_[apvId_[ihist]][stripId_[ihist]];
406  Bmin[apvId_[ihist]] = Bmin[apvId_[ihist]] < cal_->baseline_[apvId_[ihist]][stripId_[ihist]]
407  ? Bmin[apvId_[ihist]]
408  : cal_->baseline_[apvId_[ihist]][stripId_[ihist]];
409  Bmax[apvId_[ihist]] = Bmax[apvId_[ihist]] > cal_->baseline_[apvId_[ihist]][stripId_[ihist]]
410  ? Bmax[apvId_[ihist]]
411  : cal_->baseline_[apvId_[ihist]][stripId_[ihist]];
412  Bspread[apvId_[ihist]] +=
413  cal_->baseline_[apvId_[ihist]][stripId_[ihist]] * cal_->baseline_[apvId_[ihist]][stripId_[ihist]];
414  }
415 
416  // make mean values
417  for (int i = 0; i < 2; i++) {
418  if (nStrips[i] != 0) {
419  Amean[i] = Amean[i] / nStrips[i];
420  Tmean[i] = Tmean[i] / nStrips[i];
421  Rmean[i] = Rmean[i] / nStrips[i];
422  Cmean[i] = Cmean[i] / nStrips[i];
423  Omean[i] = Omean[i] / nStrips[i];
424  Mmean[i] = Mmean[i] / nStrips[i];
425  Umean[i] = Umean[i] / nStrips[i];
426  Bmean[i] = Bmean[i] / nStrips[i];
427  Smean[i] = Smean[i] / nStrips[i];
428  Kmean[i] = Kmean[i] / nStrips[i];
429 
430  Aspread[i] = Aspread[i] / nStrips[i];
431  Tspread[i] = Tspread[i] / nStrips[i];
432  Rspread[i] = Rspread[i] / nStrips[i];
433  Cspread[i] = Cspread[i] / nStrips[i];
434  Ospread[i] = Ospread[i] / nStrips[i];
435  Mspread[i] = Mspread[i] / nStrips[i];
436  Uspread[i] = Uspread[i] / nStrips[i];
437  Bspread[i] = Bspread[i] / nStrips[i];
438  Sspread[i] = Sspread[i] / nStrips[i];
439  Kspread[i] = Kspread[i] / nStrips[i];
440  }
441  }
442 
443  // fill the mean, max, min, spread, ... histograms.
444  for (int i = 0; i < 2; ++i) {
445  cal_->mean_amplitude_[i] = Amean[i];
446  cal_->mean_tail_[i] = Tmean[i];
447  cal_->mean_riseTime_[i] = Rmean[i];
448  cal_->mean_decayTime_[i] = Cmean[i];
449  cal_->mean_turnOn_[i] = Omean[i];
450  cal_->mean_peakTime_[i] = Mmean[i];
451  cal_->mean_undershoot_[i] = Umean[i];
452  cal_->mean_baseline_[i] = Bmean[i];
453  cal_->mean_smearing_[i] = Smean[i];
454  cal_->mean_chi2_[i] = Kmean[i];
455 
456  cal_->min_amplitude_[i] = Amin[i];
457  cal_->min_tail_[i] = Tmin[i];
458  cal_->min_riseTime_[i] = Rmin[i];
459  cal_->min_decayTime_[i] = Cmin[i];
460  cal_->min_turnOn_[i] = Omin[i];
461  cal_->min_peakTime_[i] = Mmin[i];
462  cal_->min_undershoot_[i] = Umin[i];
463  cal_->min_baseline_[i] = Bmin[i];
464  cal_->min_smearing_[i] = Smin[i];
465  cal_->min_chi2_[i] = Kmin[i];
466 
467  cal_->max_amplitude_[i] = Amax[i];
468  cal_->max_tail_[i] = Tmax[i];
469  cal_->max_riseTime_[i] = Rmax[i];
470  cal_->max_decayTime_[i] = Cmax[i];
471  cal_->max_turnOn_[i] = Omax[i];
472  cal_->max_peakTime_[i] = Mmax[i];
473  cal_->max_undershoot_[i] = Umax[i];
474  cal_->max_baseline_[i] = Bmax[i];
475  cal_->max_smearing_[i] = Smax[i];
476  cal_->max_chi2_[i] = Kmax[i];
477 
478  cal_->spread_amplitude_[i] = sqrt(fabs(Aspread[i] - Amean[i] * Amean[i]));
479  cal_->spread_tail_[i] = sqrt(fabs(Tspread[i] - Tmean[i] * Tmean[i]));
480  cal_->spread_riseTime_[i] = sqrt(fabs(Rspread[i] - Rmean[i] * Rmean[i]));
481  cal_->spread_decayTime_[i] = sqrt(fabs(Cspread[i] - Cmean[i] * Cmean[i]));
482  cal_->spread_turnOn_[i] = sqrt(fabs(Ospread[i] - Omean[i] * Omean[i]));
483  cal_->spread_peakTime_[i] = sqrt(fabs(Mspread[i] - Mmean[i] * Mmean[i]));
484  cal_->spread_undershoot_[i] = sqrt(fabs(Uspread[i] - Umean[i] * Umean[i]));
485  cal_->spread_baseline_[i] = sqrt(fabs(Bspread[i] - Bmean[i] * Bmean[i]));
486  cal_->spread_smearing_[i] = sqrt(fabs(Sspread[i] - Smean[i] * Smean[i]));
487  cal_->spread_chi2_[i] = sqrt(fabs(Kspread[i] - Kmean[i] * Kmean[i]));
488  }
489 
490  if (fit_function)
491  delete fit_function;
492 }
493 
494 // ------
495 void CalibrationAlgorithm::correctDistribution(TH1* histo) const {
496  // 20 events per point in the TM loop --> divide by 20 to have the amplitude of a single event readout
497  for (int iBin = 0; iBin < histo->GetNbinsX(); iBin++) {
498  histo->SetBinContent(iBin + 1, -histo->GetBinContent(iBin + 1) / 20.);
499  }
500 }
501 
502 // ----------------------------------------------------------------------------
503 float CalibrationAlgorithm::baseLine(TF1* f) {
504  float xmax = 10;
505  float baseline = 0;
506  int npoints = 0;
507  float x = f->GetXmin();
508  for (; x < xmax; x += 0.1) {
509  baseline += f->Eval(x);
510  npoints++;
511  }
512  return baseline / npoints;
513 }
514 
515 // ----------------------------------------------------------------------------
516 float CalibrationAlgorithm::turnOn(TF1* f,
517  const float& baseline) { // should happen within 100 ns in both deco and peak modes
518  float max_amplitude = f->GetMaximum();
519  float time = 10.;
520  for (; time < 100 && (f->Eval(time) - baseline) < 0.05 * (max_amplitude - baseline); time += 0.1) {
521  } // flucutation higher than 5% of the pulse height
522  return time;
523 }
524 
525 // ----------------------------------------------------------------------------
526 float CalibrationAlgorithm::decayTime(
527  TF1* f) { // if we approximate the decay to an exp(-t/tau), in one constant unit, the amplited is reduced by e^{-1}
528  float xval = f->GetMaximumX();
529  float max_amplitude = f->GetMaximum();
530  float x = xval;
531  for (; x < 1000; x = x + 0.1) {
532  if (f->Eval(x) < max_amplitude * exp(-1))
533  break;
534  }
535  return x;
536 }
CommissioningAlgorithm::anal
CommissioningAnalysis *const anal() const
Definition: CommissioningAlgorithm.h:50
sistrip::unexpectedTask_
static const char unexpectedTask_[]
Definition: ConstantsForCommissioningAnalysis.h:21
CalibrationAnalysis::minRiseTimeThresholdDeco_
static const float minRiseTimeThresholdDeco_
Definition: CalibrationAnalysis.h:119
CalibrationAnalysis::maxTurnOnThreshold_
static const float maxTurnOnThreshold_
Definition: CalibrationAnalysis.h:113
CalibrationAnalysis::maxRiseTimeThresholdDeco_
static const float maxRiseTimeThresholdDeco_
Definition: CalibrationAnalysis.h:120
CalibrationAnalysis::minPeakTimeThreshold_
static const float minPeakTimeThreshold_
Definition: CalibrationAnalysis.h:108
CalibrationAnalysis::maxDecayTimeThresholdDeco_
static const float maxDecayTimeThresholdDeco_
Definition: CalibrationAnalysis.h:116
mps_fire.i
i
Definition: mps_fire.py:429
SiStripHistoTitle
Utility class that holds histogram title.
Definition: SiStripHistoTitle.h:20
CalibrationAnalysis::maxRiseTimeThreshold_
static const float maxRiseTimeThreshold_
Definition: CalibrationAnalysis.h:111
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
CalibrationAlgorithm::correctDistribution
void correctDistribution(TH1 *) const
Definition: CalibrationAlgorithm.cc:495
CalibrationAnalysis::maxBaselineThreshold_
static const float maxBaselineThreshold_
Definition: CalibrationAnalysis.h:103
CalibrationAlgorithm::calChan_
std::vector< int > calChan_
Definition: CalibrationAlgorithm.h:43
hcalRecHitTable_cff.time
time
Definition: hcalRecHitTable_cff.py:14
sistrip
sistrip classes
Definition: EnsembleCalibrationLA.cc:9
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
me0TriggerPseudoDigis_cff.nStrips
nStrips
1.2 is to make the matching window safely the two nearest strips 0.35 is the size of an ME0 chamber i...
Definition: me0TriggerPseudoDigis_cff.py:26
CalibrationAlgorithm::cal_
CalibrationAnalysis * cal_
Definition: CalibrationAlgorithm.h:47
CalibrationAnalysis::minDecayTimeThresholdDeco_
static const float minDecayTimeThresholdDeco_
Definition: CalibrationAnalysis.h:115
CalibrationAlgorithm::extract
void extract(const std::vector< TH1 *> &) override
Definition: CalibrationAlgorithm.cc:28
CommissioningAlgorithm::extractFedKey
uint32_t extractFedKey(const TH1 *const)
Definition: CommissioningAlgorithm.cc:29
CalibrationAnalysis::minAmplitudeThreshold_
static const float minAmplitudeThreshold_
Definition: CalibrationAnalysis.h:101
CalibrationAlgorithm::histo
const Histo & histo(int &i)
Definition: CalibrationAlgorithm.h:24
CalibrationAnalysis::maxDecayTimeThreshold_
static const float maxDecayTimeThreshold_
Definition: CalibrationAnalysis.h:107
CalibrationAnalysis
Analysis for calibration runs.
Definition: CalibrationAnalysis.h:16
CalibrationAnalysis::minTurnOnThresholdDeco_
static const float minTurnOnThresholdDeco_
Definition: CalibrationAnalysis.h:121
CalibrationAnalysis::minDecayTimeThreshold_
static const float minDecayTimeThreshold_
Definition: CalibrationAnalysis.h:106
sistrip::mlCommissioning_
static const char mlCommissioning_[]
Definition: ConstantsForLogger.h:15
CalibrationAlgorithm::apvId_
std::vector< int > apvId_
Definition: CalibrationAlgorithm.h:44
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
CalibrationAnalysis::minPeakTimeThresholdDeco_
static const float minPeakTimeThresholdDeco_
Definition: CalibrationAnalysis.h:117
CommissioningAnalysis::addErrorCode
virtual void addErrorCode(const std::string &error)
Definition: CommissioningAnalysis.h:148
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
npoints
static const int npoints
Definition: NuclearInteractionFTFSimulator.h:38
CalibrationAnalysis::minBaselineThreshold_
static const float minBaselineThreshold_
Definition: CalibrationAnalysis.h:102
HLT_2024v12_cff.Tmin
Tmin
Definition: HLT_2024v12_cff.py:18319
CommissioningAnalysis::fedKey
const uint32_t & fedKey() const
Definition: CommissioningAnalysis.h:134
CalibrationAlgorithm::histo_
std::vector< Histo > histo_
Definition: CalibrationAlgorithm.h:41
CalibrationAnalysis::maxTurnOnThresholdDeco_
static const float maxTurnOnThresholdDeco_
Definition: CalibrationAnalysis.h:122
CalibrationAlgorithm::stripId_
std::vector< int > stripId_
Definition: CalibrationAlgorithm.h:42
CalibrationAnalysis::minRiseTimeThreshold_
static const float minRiseTimeThreshold_
Definition: CalibrationAnalysis.h:110
TrackerOfflineValidation_Dqm_cff.xmax
xmax
Definition: TrackerOfflineValidation_Dqm_cff.py:11
fdeconv
double fdeconv(double *x, double *par)
Definition: SiStripPulseShape.cc:101
CommissioningAlgorithm::Histo
std::pair< TH1 *, std::string > Histo
Definition: CommissioningAlgorithm.h:25
CalibrationAnalysis::minTurnOnThreshold_
static const float minTurnOnThreshold_
Definition: CalibrationAnalysis.h:112
CalibrationAnalysis::maxChi2Threshold_
static const float maxChi2Threshold_
Definition: CalibrationAnalysis.h:104
combine.histos
histos
Definition: combine.py:4
CalibrationAnalysis::maxPeakTimeThresholdDeco_
static const float maxPeakTimeThresholdDeco_
Definition: CalibrationAnalysis.h:118
fpeak
double fpeak(double *x, double *par)
Definition: SiStripPulseShape.cc:66
CalibrationAlgorithm::turnOn
float turnOn(TF1 *, const float &)
Definition: CalibrationAlgorithm.cc:516
CalibrationAnalysis::maxPeakTimeThreshold_
static const float maxPeakTimeThreshold_
Definition: CalibrationAnalysis.h:109
Gflash::Rmax
const double Rmax[kNumberCalorimeter]
Definition: GflashNameSpace.h:31
CommissioningAnalysis
Abstract base for derived classes that provide analysis of commissioning histograms.
Definition: CommissioningAnalysis.h:18
conversion_template_cfg.anal
anal
Definition: conversion_template_cfg.py:16
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:128
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
Gflash::Rmin
const double Rmin[kNumberCalorimeter]
Definition: GflashNameSpace.h:30
Generated for CMSSW Reference Manual by doxygen 1.8.14