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 <iostream>
13 #include <sstream>
14 #include <iomanip>
15 #include <cmath>
16 
17 using namespace sistrip;
18 
19 // ----------------------------------------------------------------------------
20 //
21 CalibrationAlgorithm::CalibrationAlgorithm( const edm::ParameterSet & pset, CalibrationAnalysis* const anal )
22  : CommissioningAlgorithm(anal),
23  deconv_fitter_(nullptr),
24  peak_fitter_(nullptr),
25  cal_(nullptr)
26 {
27 
28  deconv_fitter_ = new TF1("deconv_fitter",fdeconv,0,200,6);
29  deconv_fitter_->SetParLimits(0,1,70); //x
30  deconv_fitter_->SetParLimits(1,0,40); //z = tau
31  deconv_fitter_->SetParLimits(4,15,55); // turn-on-time t_0
32  deconv_fitter_->SetParameters(2,25,0.8000,50,40);
33 
34  peak_fitter_ = new TF1("peak_fitter",fpeak,0,200,5);
35  peak_fitter_->SetParLimits(0,1,70); //x
36  peak_fitter_->SetParLimits(1,20,70);//z = tau
37  peak_fitter_->SetParLimits(4,15,35); //turn-on-time t_0
38  peak_fitter_->SetParameters(17,50,0,5000,20);
39 }
40 
41 // ----------------------------------------------------------------------------
42 //
43 void CalibrationAlgorithm::extract( const std::vector<TH1*>& histos) {
44  if ( !anal() ) {
45  edm::LogWarning(mlCommissioning_)
46  << "[CalibrationAlgorithm::" << __func__ << "]"
47  << " NULL pointer to base Analysis object!";
48  return;
49  }
50 
51  CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>( anal() );
52  cal_ = dynamic_cast<CalibrationAnalysis*>( tmp );
53  if ( !cal_ ) {
54  edm::LogWarning(mlCommissioning_)
55  << "[CalibrationAlgorithm::" << __func__ << "]"
56  << " NULL pointer to derived Analysis object!";
57  return;
58  }
59 
60  // Check number of histograms
61  if ( histos.size() != 32 && histos.size() !=2 ) {
62  cal_->addErrorCode(sistrip::numberOfHistos_);
63  }
64 
65  // Extract FED key from histo title
66  if ( !histos.empty() ) { cal_->fedKey( extractFedKey( histos.front() ) ); }
67 
68  // Extract histograms
69  std::vector<TH1*>::const_iterator ihis = histos.begin();
70  unsigned int cnt = 0;
71  for ( ; ihis != histos.end(); ihis++,cnt++ ) {
72 
73  // Check for NULL pointer
74  if ( !(*ihis) ) { continue; }
75 
76  // Check name
77  SiStripHistoTitle title( (*ihis)->GetName() );
78  if ( title.runType() != sistrip::CALIBRATION &&
79  title.runType() != sistrip::CALIBRATION_DECO &&
80  title.runType() != sistrip::CALIBRATION_SCAN &&
81  title.runType() != sistrip::CALIBRATION_SCAN_DECO ) {
82  cal_->addErrorCode(sistrip::unexpectedTask_);
83  continue;
84  }
85 
86  cal_->isScan_ = ( title.runType() == sistrip::CALIBRATION_SCAN ||
87  title.runType() == sistrip::CALIBRATION_SCAN_DECO );
88 
89  // Extract calibration histo
90  histo_[cnt].first = *ihis;
91  histo_[cnt].second = (*ihis)->GetTitle();
92  }
93 
94 }
95 
96 // ----------------------------------------------------------------------------
97 //
98 void CalibrationAlgorithm::analyse() {
99 
100  if ( !cal_ ) {
101  edm::LogWarning(mlCommissioning_)
102  << "[CalibrationAlgorithm::" << __func__ << "]"
103  << " NULL pointer to derived Analysis object!";
104  return;
105  }
106 
107  float Amean[2] = {0.,0.};
108  float Amin[2] = {2000.,2000.};
109  float Amax[2] = {0.,0.};
110  float Aspread[2] = {0.,0.};
111  float Tmean[2] = {0.,0.};
112  float Tmin[2] = {2000.,2000.};
113  float Tmax[2] = {0.,0.};
114  float Tspread[2] = {0.,0.};
115  float Rmean[2] = {0.,0.};
116  float Rmin[2] = {2000.,2000.};
117  float Rmax[2] = {0.,0.};
118  float Rspread[2] = {0.,0.};
119  float Cmean[2] = {0.,0.};
120  float Cmin[2] = {2000.,2000.};
121  float Cmax[2] = {0.,0.};
122  float Cspread[2] = {0.,0.};
123  float Smean[2] = {0.,0.};
124  float Smin[2] = {2000.,2000.};
125  float Smax[2] = {0.,0.};
126  float Sspread[2] = {0.,0.};
127  float Kmean[2] = {0.,0.};
128  float Kmin[2] = {2000000.,2000000.};
129  float Kmax[2] = {0.,0.};
130  float Kspread[2] = {0.,0.};
131  // turnOn
132  float Omean[2] = {0.,0.};
133  float Omin[2] = {2000.,2000.};
134  float Omax[2] = {0.,0.};
135  float Ospread[2] = {0.,0.};
136  // maximum
137  float Mmean[2] = {0.,0.};
138  float Mmin[2] = {2000.,2000.};
139  float Mmax[2] = {0.,0.};
140  float Mspread[2] = {0.,0.};
141  // undershoot
142  float Umean[2] = {0.,0.};
143  float Umin[2] = {2000.,2000.};
144  float Umax[2] = {0.,0.};
145  float Uspread[2] = {0.,0.};
146  // baseline
147  float Bmean[2] = {0.,0.};
148  float Bmin[2] = {2000.,2000.};
149  float Bmax[2] = {0.,0.};
150  float Bspread[2] = {0.,0.};
151 
152  unsigned int upperLimit = cal_->isScan_ ? 2 : 32;
153  float nStrips = cal_->isScan_ ? 1. : 16.;
154  for(unsigned int i=0;i<upperLimit;++i) {
155  if ( !histo_[i].first ) {
156  edm::LogWarning(mlCommissioning_)
157  << " NULL pointer to histogram " << i << "!";
158  return;
159  }
160 
161  // determine which APV and strip is being looked at.
162  std::string title = histo_[i].first->GetName();
163  int apv = 0;
164  int strip = 0;
165  if(title.find("STRIP")!=std::string::npos && title.find("Apv")!=std::string::npos) {
166  strip = atoi(title.c_str()+title.find("STRIP")+6);
167  apv = (atoi(title.c_str()+title.find("Apv")+3))%2;
168  } else {
169  strip = (atoi(title.c_str()+title.find_last_of("_")+1))%16;
170  apv = (atoi(title.c_str()+title.find_last_of("_")+1))/16;
171  if(title.find("Apv")!=std::string::npos) {
172  apv = (atoi(title.c_str()+title.find("Apv")+3))%2;
173  strip = strip*8 + cal_->calchan_;
174  } else {
175  edm::LogWarning(mlCommissioning_)
176  << " Malformed histogram title! Strip/APV not retreived: "
177  << title;
178  }
179  }
180 
181  // rescale the plot
182  correctDistribution(histo_[i].first);
183 
185  TF1* fit = fitPulse(histo_[i].first);
186  float maximum_ampl = fit->GetMaximum();
187  float baseline = fit->Eval(10);
188  cal_->amplitude_[apv][strip] = maximum_ampl - baseline;
189 
190  // rise time
191  float peak_time = fit->GetMaximumX();
192  float turn_on_time = turnOn(fit);
193  float rise_time = peak_time - turn_on_time;
194  cal_->riseTime_[apv][strip] = rise_time;
195 
196  //turn-on
197  cal_->turnOn_[apv][strip] = turn_on_time;
198 
199  //maximum
200  cal_->maximum_[apv][strip] = peak_time;
201 
202  //undershoot
203  if (cal_->deconv_)
204  cal_->undershoot_[apv][strip] = 100*(fit->GetMinimum()-baseline)/(maximum_ampl - baseline);
205  else
206  cal_->undershoot_[apv][strip] = 0;
207 
208  //baseline
209  cal_->baseline_[apv][strip] = 100 * baseline/(maximum_ampl - baseline);
210 
211 
212  // tail 125 ns after the maximum
213  int lastBin = histo_[i].first->FindBin(peak_time + 125);
214  if(lastBin>histo_[i].first->GetNbinsX()-4) lastBin = histo_[i].first->GetNbinsX()-4;
215  if(histo_[i].first->GetMaximum()!=0)
216  cal_->tail_[apv][strip] = 100*(histo_[i].first->GetBinContent(lastBin)-baseline) / (maximum_ampl - baseline);
217  else
218  cal_->tail_[apv][strip] = 100;
219 
220  cal_->timeConstant_[apv][strip] = fit->GetParameter(1);
221  cal_->smearing_[apv][strip] = 0;
222  cal_->chi2_[apv][strip] = fit->GetChisquare();
223 
224 
225  //compute mean, max, min, spread
226  Amean[apv] += cal_->amplitude_[apv][strip]/nStrips;
227  Amin[apv] = Amin[apv]<cal_->amplitude_[apv][strip] ? Amin[apv] : cal_->amplitude_[apv][strip];
228  Amax[apv] = Amax[apv]>cal_->amplitude_[apv][strip] ? Amax[apv] : cal_->amplitude_[apv][strip];
229  Aspread[apv] += cal_->amplitude_[apv][strip]*cal_->amplitude_[apv][strip]/nStrips;
230 
231  Tmean[apv] += cal_->tail_[apv][strip]/nStrips;
232  Tmin[apv] = Tmin[apv]<cal_->tail_[apv][strip] ? Tmin[apv] : cal_->tail_[apv][strip];
233  Tmax[apv] = Tmax[apv]>cal_->tail_[apv][strip] ? Tmax[apv] : cal_->tail_[apv][strip];
234  Tspread[apv] += cal_->tail_[apv][strip]*cal_->tail_[apv][strip]/nStrips;
235 
236  Rmean[apv] += cal_->riseTime_[apv][strip]/nStrips;
237  Rmin[apv] = Rmin[apv]<cal_->riseTime_[apv][strip] ? Rmin[apv] : cal_->riseTime_[apv][strip];
238  Rmax[apv] = Rmax[apv]>cal_->riseTime_[apv][strip] ? Rmax[apv] : cal_->riseTime_[apv][strip];
239  Rspread[apv] += cal_->riseTime_[apv][strip]*cal_->riseTime_[apv][strip]/nStrips;
240 
241  Cmean[apv] += cal_->timeConstant_[apv][strip]/nStrips;
242  Cmin[apv] = Cmin[apv]<cal_->timeConstant_[apv][strip] ? Cmin[apv] : cal_->timeConstant_[apv][strip];
243  Cmax[apv] = Cmax[apv]>cal_->timeConstant_[apv][strip] ? Cmax[apv] : cal_->timeConstant_[apv][strip];
244  Cspread[apv] += cal_->timeConstant_[apv][strip]*cal_->timeConstant_[apv][strip]/nStrips;
245 
246  Smean[apv] += cal_->smearing_[apv][strip]/nStrips;
247  Smin[apv] = Smin[apv]<cal_->smearing_[apv][strip] ? Smin[apv] : cal_->smearing_[apv][strip];
248  Smax[apv] = Smax[apv]>cal_->smearing_[apv][strip] ? Smax[apv] : cal_->smearing_[apv][strip];
249  Sspread[apv] += cal_->smearing_[apv][strip]*cal_->smearing_[apv][strip]/nStrips;
250 
251  Kmean[apv] += cal_->chi2_[apv][strip]/nStrips;
252  Kmin[apv] = Kmin[apv]<cal_->chi2_[apv][strip] ? Kmin[apv] : cal_->chi2_[apv][strip];
253  Kmax[apv] = Kmax[apv]>cal_->chi2_[apv][strip] ? Kmax[apv] : cal_->chi2_[apv][strip];
254  Kspread[apv] += cal_->chi2_[apv][strip]*cal_->chi2_[apv][strip]/nStrips;
255 
256  Omean[apv] += cal_->turnOn_[apv][strip]/nStrips;
257  Omin[apv] = Omin[apv]<cal_->turnOn_[apv][strip] ? Omin[apv] : cal_->turnOn_[apv][strip];
258  Omax[apv] = Omax[apv]>cal_->turnOn_[apv][strip] ? Omax[apv] : cal_->turnOn_[apv][strip];
259  Ospread[apv] += cal_->turnOn_[apv][strip]*cal_->turnOn_[apv][strip]/nStrips;
260 
261  Mmean[apv] += cal_->maximum_[apv][strip]/nStrips;
262  Mmin[apv] = Mmin[apv]<cal_->maximum_[apv][strip] ? Mmin[apv] : cal_->maximum_[apv][strip];
263  Mmax[apv] = Mmax[apv]>cal_->maximum_[apv][strip] ? Mmax[apv] : cal_->maximum_[apv][strip];
264  Mspread[apv] += cal_->maximum_[apv][strip]*cal_->maximum_[apv][strip]/nStrips;
265 
266  Umean[apv] += cal_->undershoot_[apv][strip]/nStrips;
267  Umin[apv] = Umin[apv]<cal_->undershoot_[apv][strip] ? Umin[apv] : cal_->undershoot_[apv][strip];
268  Umax[apv] = Umax[apv]>cal_->undershoot_[apv][strip] ? Umax[apv] : cal_->undershoot_[apv][strip];
269  Uspread[apv] += cal_->undershoot_[apv][strip]*cal_->undershoot_[apv][strip]/nStrips;
270 
271  Bmean[apv] += cal_->baseline_[apv][strip]/nStrips;
272  Bmin[apv] = Bmin[apv]<cal_->baseline_[apv][strip] ? Bmin[apv] : cal_->baseline_[apv][strip];
273  Bmax[apv] = Bmax[apv]>cal_->baseline_[apv][strip] ? Bmax[apv] : cal_->baseline_[apv][strip];
274  Bspread[apv] += cal_->baseline_[apv][strip]*cal_->baseline_[apv][strip]/nStrips;
275  }
276 
277  // fill the mean, max, min, spread, ... histograms.
278  for(int i=0;i<2;++i) {
279  cal_->mean_amplitude_[i] = Amean[i];
280  cal_->mean_tail_[i] = Tmean[i];
281  cal_->mean_riseTime_[i] = Rmean[i];
282  cal_->mean_timeConstant_[i] = Cmean[i];
283  cal_->mean_turnOn_[i] = Omean[i];
284  cal_->mean_maximum_[i] = Mmean[i];
285  cal_->mean_undershoot_[i] = Umean[i];
286  cal_->mean_baseline_[i] = Bmean[i];
287  cal_->mean_smearing_[i] = Smean[i];
288  cal_->mean_chi2_[i] = Kmean[i];
289  cal_->min_amplitude_[i] = Amin[i];
290  cal_->min_tail_[i] = Tmin[i];
291  cal_->min_riseTime_[i] = Rmin[i];
292  cal_->min_timeConstant_[i] = Cmin[i];
293  cal_->min_turnOn_[i] = Omin[i];
294  cal_->min_maximum_[i] = Mmin[i];
295  cal_->min_undershoot_[i] = Umin[i];
296  cal_->min_baseline_[i] = Bmin[i];
297  cal_->min_smearing_[i] = Smin[i];
298  cal_->min_chi2_[i] = Kmin[i];
299  cal_->max_amplitude_[i] = Amax[i];
300  cal_->max_tail_[i] = Tmax[i];
301  cal_->max_riseTime_[i] = Rmax[i];
302  cal_->max_timeConstant_[i] = Cmax[i];
303  cal_->max_turnOn_[i] = Omax[i];
304  cal_->max_maximum_[i] = Mmax[i];
305  cal_->max_undershoot_[i] = Umax[i];
306  cal_->max_baseline_[i] = Bmax[i];
307  cal_->max_smearing_[i] = Smax[i];
308  cal_->max_chi2_[i] = Kmax[i];
309  cal_->spread_amplitude_[i] = sqrt(fabs(Aspread[i]-Amean[i]*Amean[i]));
310  cal_->spread_tail_[i] = sqrt(fabs(Tspread[i]-Tmean[i]*Tmean[i]));
311  cal_->spread_riseTime_[i] = sqrt(fabs(Rspread[i]-Rmean[i]*Rmean[i]));
312  cal_->spread_timeConstant_[i] = sqrt(fabs(Cspread[i]-Cmean[i]*Cmean[i]));
313  cal_->spread_turnOn_[i] = sqrt(fabs(Ospread[i]-Omean[i]*Omean[i]));
314  cal_->spread_maximum_[i] = sqrt(fabs(Mspread[i]-Mmean[i]*Mmean[i]));
315  cal_->spread_undershoot_[i] = sqrt(fabs(Uspread[i]-Umean[i]*Umean[i]));
316  cal_->spread_baseline_[i] = sqrt(fabs(Bspread[i]-Bmean[i]*Bmean[i]));
317  cal_->spread_smearing_[i] = sqrt(fabs(Sspread[i]-Smean[i]*Smean[i]));
318  cal_->spread_chi2_[i] = sqrt(fabs(Kspread[i]-Kmean[i]*Kmean[i]));
319  }
320 }
321 
322 // ----------------------------------------------------------------------------
323 //
324 void CalibrationAlgorithm::correctDistribution( TH1* histo ) const
325 {
326  // return the curve
327  histo->Scale(-1);
328  if ( cal_ ) {
329  if( cal_->isScan_ ){
330  histo->Scale(1/16.);
331  edm::LogWarning(mlCommissioning_)
332  << "CalibrationAlgorithm::correctDistribution: isScan_ == true!!! ";
333  }
334  }
335 }
336 
337 // ----------------------------------------------------------------------------
338 //
339 TF1* CalibrationAlgorithm::fitPulse( TH1* histo,
340  float rangeLow,
341  float rangeHigh )
342 {
343  if(!cal_) return nullptr;
344  if(!histo) return nullptr;
345  float noise = 4.;
346  float N = round(histo->GetMaximum()/125.);
347  float error = sqrt(2*N)*noise;
348 
349  for(int i=1;i<=histo->GetNbinsX();++i) {
350  histo->SetBinError(i,error);
351  }
352 
353  // let's help our fit with a bit of parameter guessing
354  float maximum = histo->GetMaximum();
355  float turn_on = histo->GetBinCenter((histo->FindFirstBinAbove(0.1 * maximum) - 2));
356 
357  if (cal_->deconv_) {
358 
359  deconv_fitter_->SetParameter(4, turn_on);
360 
361  if(rangeLow>rangeHigh)
362  histo->Fit(deconv_fitter_,"QMR");
363  else
364  histo->Fit(deconv_fitter_,"0QMR","",rangeLow,rangeHigh);
365 
366  return deconv_fitter_;
367 
368  }
369  else {
370  peak_fitter_->SetParameter(4, turn_on);
371  if(rangeLow>rangeHigh)
372  histo->Fit(peak_fitter_,"MRQ");
373  else
374  histo->Fit(peak_fitter_,"0QMR","",rangeLow,rangeHigh);
375  return peak_fitter_;
376  }
377 }
378 
379 // ----------------------------------------------------------------------------
380 //
381 float CalibrationAlgorithm::maximum( TH1* h ) {
382  int bin = h->GetMaximumBin();
383  // fit around the maximum with the detector response and take the max from the fit
384  TF1* fit = fitPulse(h,h->GetBinCenter(bin)-25,h->GetBinCenter(bin)+25);
385  return fit->GetMaximumX();
386 }
387 
388 // ----------------------------------------------------------------------------
389 //
390 float CalibrationAlgorithm::turnOn(TF1* f){
391  float max_amplitude = f->GetMaximum();
392  float time = 5.;
393  float base = f->GetMinimum(0,50);
394  for( ; time < 50 && (f->Eval(time) - base) < 0.01 * (max_amplitude - base); time += 0.1) {}
395  return time - 0.05;
396 }
CalibrationAlgorithm::fitPulse
TF1 * fitPulse(TH1 *, float rangeLow=0, float rangeHigh=-1)
Definition: CalibrationAlgorithm.cc:339
sistrip::unexpectedTask_
static const char unexpectedTask_[]
Definition: ConstantsForCommissioningAnalysis.h:21
mps_fire.i
i
Definition: mps_fire.py:269
CommissioningAnalysis::fedKey
const uint32_t & fedKey() const
Definition: CommissioningAnalysis.h:136
h
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
Definition: L1TUtmAlgorithmRcd.h:4
CalibrationAlgorithm::histo
const Histo & histo(int i) const
Definition: CalibrationAlgorithm.h:26
SiStripHistoTitle
Utility class that holds histogram title.
Definition: SiStripHistoTitle.h:20
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
sistrip::numberOfHistos_
static const char numberOfHistos_[]
Definition: ConstantsForCommissioningAnalysis.h:16
CalibrationAlgorithm::extract
void extract(const std::vector< TH1 * > &) override
Definition: CalibrationAlgorithm.cc:43
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
nullptr
#define nullptr
sistrip
sistrip classes
Definition: EnsembleCalibrationLA.cc:12
CalibrationAlgorithm::cal_
CalibrationAnalysis * cal_
Definition: CalibrationAlgorithm.h:56
CalibrationAnalysis
Analysis for calibration runs.
Definition: CalibrationAnalysis.h:16
sistrip::mlCommissioning_
static const char mlCommissioning_[]
Definition: ConstantsForLogger.h:15
OfflinePixel3DPrimaryVertices_cfi.Tmin
Tmin
Definition: OfflinePixel3DPrimaryVertices_cfi.py:25
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
CommissioningAlgorithm::extractFedKey
uint32_t extractFedKey(const TH1 *const )
Definition: CommissioningAlgorithm.cc:31
CommissioningAnalysis::addErrorCode
virtual void addErrorCode(const std::string &error)
Definition: CommissioningAnalysis.h:150
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
runEdmFileComparison.base
base
Make Sure CMSSW is Setup ##.
Definition: runEdmFileComparison.py:159
stringResolutionProvider_cfi.bin
bin
set the eta bin as selection string.
Definition: stringResolutionProvider_cfi.py:11
N
#define N
Definition: blowfish.cc:9
fdeconv
double fdeconv(double *x, double *par)
Definition: SiStripPulseShape.cc:102
CalibrationAlgorithm::correctDistribution
void correctDistribution(TH1 *) const
Definition: CalibrationAlgorithm.cc:324
digitizers_cfi.strip
strip
Definition: digitizers_cfi.py:19
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
fftjetcommon_cfi.title
title
Definition: fftjetcommon_cfi.py:32
fpeak
double fpeak(double *x, double *par)
Definition: SiStripPulseShape.cc:87
Gflash::Rmax
const double Rmax[kNumberCalorimeter]
Definition: GflashNameSpace.h:30
CommissioningAnalysis
Abstract base for derived classes that provide analysis of commissioning histograms.
Definition: CommissioningAnalysis.h:18
Gflash::Rmin
const double Rmin[kNumberCalorimeter]
Definition: GflashNameSpace.h:29
CommissioningAlgorithm::anal
CommissioningAnalysis *const anal() const
Definition: CommissioningAlgorithm.h:55
Generated for CMSSW Reference Manual by doxygen 1.8.11