CMS 3D CMS Logo

SamplingAlgorithm.cc
Go to the documentation of this file.
7 #include "TProfile.h"
8 #include "TF1.h"
9 #include <iostream>
10 #include <sstream>
11 #include <iomanip>
12 #include <cmath>
13 
14 using namespace sistrip;
15 
16 // ----------------------------------------------------------------------------
17 //
20  histo_(nullptr, ""),
21  deconv_fitter_(nullptr),
22  peak_fitterA_(nullptr),
23  peak_fitterB_(nullptr),
24  latencyCode_(latencyCode),
25  samp_(nullptr) {
26  peak_fitterA_ = new TF1("peak_fitterA", fpeak_convoluted, -4800, 0, 5);
27  peak_fitterA_->SetNpx(2000);
28  peak_fitterA_->FixParameter(0, 0);
29  peak_fitterA_->SetParLimits(1, 0, 4800);
30  peak_fitterA_->SetParLimits(2, 0, 20);
31  peak_fitterA_->FixParameter(3, 50);
32  peak_fitterA_->SetParLimits(4, 0, 75);
33  peak_fitterA_->SetParameters(0., 1250, 10, 50, 10);
34 
35  peak_fitterB_ = new TF1("peak_fitterB", fpeak_convoluted, -100, 100, 5);
36  peak_fitterB_->SetNpx(200);
37  peak_fitterB_->FixParameter(0, 0);
38  peak_fitterB_->SetParLimits(1, -100, 100);
39  peak_fitterB_->SetParLimits(2, 0, 20);
40  peak_fitterB_->FixParameter(3, 50);
41  peak_fitterB_->SetParLimits(4, 0, 75);
42  peak_fitterB_->SetParameters(0., -50, 10, 50, 10);
43 
44  deconv_fitter_ = new TF1("deconv_fitter", fdeconv_convoluted, -50, 50, 5);
45  deconv_fitter_->SetNpx(1000);
46  deconv_fitter_->FixParameter(0, 0);
47  deconv_fitter_->SetParLimits(1, -50, 50);
48  deconv_fitter_->SetParLimits(2, 0, 200);
49  deconv_fitter_->SetParLimits(3, 5, 100);
50  deconv_fitter_->FixParameter(3, 50);
51  deconv_fitter_->SetParLimits(4, 0, 75);
52  deconv_fitter_->SetParameters(0., -2.82, 0.96, 50, 20);
53 }
54 
55 // ----------------------------------------------------------------------------
56 //
57 void SamplingAlgorithm::extract(const std::vector<TH1*>& histos) {
58  if (!anal()) {
59  edm::LogWarning(mlCommissioning_) << "[SamplingAlgorithm::" << __func__ << "]"
60  << " NULL pointer to Analysis object!";
61  return;
62  }
63 
65  samp_ = dynamic_cast<SamplingAnalysis*>(tmp);
66  if (!samp_) {
67  edm::LogWarning(mlCommissioning_) << "[SamplingAlgorithm::" << __func__ << "]"
68  << " NULL pointer to derived Analysis object!";
69  return;
70  }
71 
72  // Check
73  if (histos.size() != 1 && histos.size() != 2) {
75  }
76 
77  // Extract FED key from histo title
78  if (!histos.empty()) {
79  samp_->fedKey(extractFedKey(histos.front()));
80  }
81 
82  // Extract
83  std::vector<TH1*>::const_iterator ihis = histos.begin();
84  for (; ihis != histos.end(); ihis++) {
85  // Check pointer
86  if (!(*ihis)) {
87  edm::LogWarning(mlCommissioning_) << " NULL pointer to histogram!";
88  continue;
89  }
90 
91  // Check name
92  SiStripHistoTitle title((*ihis)->GetName());
93  if (title.runType() != sistrip::APV_LATENCY && title.runType() != sistrip::FINE_DELAY) {
95  continue;
96  }
97  // Set the mode for later fits
98  samp_->runType_ = title.runType();
99 
100  // Set the granularity
101  samp_->granularity_ = title.granularity();
102 
103  // Extract timing histo
104  if (title.extraInfo().find(sistrip::extrainfo::clusterCharge_) != std::string::npos) {
105  histo_.first = *ihis;
106  histo_.second = (*ihis)->GetName();
107  }
108  }
109 }
110 
111 // ----------------------------------------------------------------------------
112 //
114  if (!samp_) {
115  edm::LogWarning(mlCommissioning_) << "[SamplingAlgorithm::" << __func__ << "]"
116  << " NULL pointer to derived Analysis object!";
117  return;
118  }
119 
120  TProfile* prof = (TProfile*)(histo_.first);
121  if (!prof) {
122  edm::LogWarning(mlCommissioning_) << " NULL pointer to histogram!";
123  return;
124  }
125 
126  // set the right error mode: rms
127  prof->SetErrorOption(" ");
128 
129  //that should not be needed, but it seems histos are stored with error option " " and errors "s" in all cases.
130  //it MUST be removed if the DQM (?) bug is solved
131  for (int i = 0; i < prof->GetNbinsX(); ++i) {
132  if (prof->GetBinEntries(i) > 0)
133  prof->SetBinError(i, prof->GetBinError(i) / sqrt(prof->GetBinEntries(i)));
134  }
135 
136  // prune the profile
137  pruneProfile(prof);
138 
139  // correct for the binning
140  correctBinning(prof);
141 
142  // correct for clustering effects
143  correctProfile(prof, samp_->sOnCut_);
144 
145  // fit depending on the mode
147  // initialize the fit (overal latency)
148  float max = prof->GetBinCenter(prof->GetMaximumBin());
149  float ampl = prof->GetMaximum();
150  peak_fitterA_->SetParameters(0., 50 - max, ampl / 20., 50, 10);
151 
152  // fit
153  if (prof->Fit(peak_fitterA_, "Q") == 0)
154  prof->Fit(peak_fitterA_, "QEM");
155 
156  // Set monitorables
157  samp_->max_ = peak_fitterA_->GetMaximumX();
158  samp_->error_ = peak_fitterA_->GetParError(1);
159 
160  } else { // sistrip::FINE_DELAY
161 
162  // initialize the fit (overal latency)
163  float max = prof->GetBinCenter(prof->GetMaximumBin());
164  float ampl = prof->GetMaximum();
165  deconv_fitter_->SetParameters(0., -max, ampl / 10., 50, 20);
166  peak_fitterB_->SetParameters(0., 50 - max, ampl / 20., 50, 10);
167  if (latencyCode_ & 0x80) { // deconv mode
168  // fit
169  if (prof->Fit(deconv_fitter_, "Q") == 0)
170  prof->Fit(deconv_fitter_, "QEM");
171  // Set monitorables
172  samp_->max_ = deconv_fitter_->GetMaximumX();
173  samp_->error_ = deconv_fitter_->GetParError(1);
174  } else { // peak mode
175  // fit
176  if (prof->Fit(peak_fitterB_, "Q") == 0)
177  prof->Fit(peak_fitterB_, "QEM");
178  // Set monitorables
179  samp_->max_ = peak_fitterB_->GetMaximumX();
180  samp_->error_ = peak_fitterB_->GetParError(1);
181  }
182  }
183 }
184 
185 // ----------------------------------------------------------------------------
186 //
188  // loop over bins to find the max stat
189  uint32_t nbins = profile->GetNbinsX();
190  uint32_t max = 0;
191  for (uint32_t bin = 1; bin <= nbins; ++bin)
192  max = max < profile->GetBinEntries(bin) ? uint32_t(profile->GetBinEntries(bin)) : max;
193  // loop over bins to clean
194  uint32_t min = max / 10;
195  for (uint32_t bin = 1; bin <= nbins; ++bin)
196  if (profile->GetBinEntries(bin) < min) {
197  profile->SetBinContent(bin, 0.);
198  profile->SetBinError(bin, 0.);
199  }
200 }
201 
202 // ----------------------------------------------------------------------------
203 //
204 void SamplingAlgorithm::correctBinning(TProfile* prof) const {
205  prof->GetXaxis()->SetLimits(prof->GetXaxis()->GetXmin() - prof->GetBinWidth(1) / 2.,
206  prof->GetXaxis()->GetXmax() - prof->GetBinWidth(1) / 2.);
207 }
208 
209 // ----------------------------------------------------------------------------
210 //
211 void SamplingAlgorithm::correctProfile(TProfile* profile, float SoNcut) const {
212  if (!samp_) {
213  return;
214  }
215  uint32_t nbins = profile->GetNbinsX();
216  float min = samp_->limit(SoNcut);
217  for (uint32_t bin = 1; bin <= nbins; ++bin)
218  if (profile->GetBinContent(bin) < min) {
219  profile->SetBinContent(bin, 0.);
220  profile->SetBinError(bin, 0.);
221  profile->SetBinEntries(bin, 0);
222  } else {
223  profile->SetBinContent(
224  bin, profile->GetBinEntries(bin) * samp_->correctMeasurement(profile->GetBinContent(bin), SoNcut));
225  }
226 }
void correctProfile(TProfile *profile, float SoNcut=3.) const
CommissioningAnalysis *const anal() const
static const char unexpectedTask_[]
Analysis for latency run.
double fdeconv_convoluted(double *x, double *par)
void extract(const std::vector< TH1 *> &) override
Utility class that holds histogram title.
sistrip::RunType runType_
static const char numberOfHistos_[]
void pruneProfile(TProfile *profile) const
void correctBinning(TProfile *prof) const
sistrip classes
SamplingAnalysis * samp_
double fpeak_convoluted(double *x, double *par)
uint32_t extractFedKey(const TH1 *const)
static const char mlCommissioning_[]
T sqrt(T t)
Definition: SSEVec.h:19
virtual void addErrorCode(const std::string &error)
const uint32_t & fedKey() const
float limit(float SoNcut) const
static const char clusterCharge_[]
float correctMeasurement(float mean, float SoNcut=3.) const
void analyse() override
histos
Definition: combine.py:4
Abstract base for derived classes that provide analysis of commissioning histograms.
Log< level::Warning, false > LogWarning
tmp
align.sh
Definition: createJobs.py:716
sistrip::Granularity granularity_