CMS 3D CMS Logo
PedsFullNoiseAlgorithm.cc
Go to the documentation of this file.
1 #include "DQM/SiStripCommissioningAnalysis/interface/PedsFullNoiseAlgorithm.h"
2 #include "CondFormats/SiStripObjects/interface/PedsFullNoiseAnalysis.h"
3 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
4 #include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h"
5 #include "FWCore/MessageLogger/interface/MessageLogger.h"
6 
7 #include "TProfile.h"
8 #include "TH1.h"
9 #include "TH2.h"
10 #include "TF1.h"
11 #include "TFitResult.h"
12 #include "TMath.h"
13 #include "Math/DistFunc.h"
14 #include "Math/ProbFuncMathCore.h"
15 #include "Fit/BinData.h"
16 #include "HFitInterface.h"
17 #include "Math/GoFTest.h"
18 #include <iostream>
19 #include <iomanip>
20 #include <cmath>
21 
22 using namespace sistrip;
23 using namespace std;
24 // ----------------------------------------------------------------------------
25 //
26 PedsFullNoiseAlgorithm::PedsFullNoiseAlgorithm( const edm::ParameterSet & pset, PedsFullNoiseAnalysis* const anal ):
27  CommissioningAlgorithm(anal),
28  hPeds_(nullptr,""),
29  hNoise_(nullptr,""),
30  hNoise2D_(nullptr,""),
31  maxDriftResidualCut_(pset.getParameter<double>("MaxDriftResidualCut")),
32  minStripNoiseCut_(pset.getParameter<double>("MinStripNoiseCut")),
33  maxStripNoiseCut_(pset.getParameter<double>("MaxStripNoiseCut")),
34  maxStripNoiseSignificanceCut_(pset.getParameter<double>("MaxStripNoiseSignificanceCut")),
35  adProbabCut_(pset.getParameter<double>("AdProbabCut")),
36  ksProbabCut_(pset.getParameter<double>("KsProbabCut")),
37  generateRandomHisto_(pset.getParameter<bool>("GenerateRandomHisto")),
38  jbProbabCut_(pset.getParameter<double>("JbProbabCut")),
39  chi2ProbabCut_(pset.getParameter<double>("Chi2ProbabCut")),
40  kurtosisCut_(pset.getParameter<double>("KurtosisCut")),
41  integralTailCut_(pset.getParameter<double>("IntegralTailCut")),
42  integralNsigma_(pset.getParameter<int>("IntegralNsigma")),
43  ashmanDistance_(pset.getParameter<double>("AshmanDistance")),
44  amplitudeRatio_(pset.getParameter<double>("AmplitudeRatio"))
45 {
46  LogDebug(mlCommissioning_)
47  << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
48  << " Set maximum drift of the mean value to: " << maxDriftResidualCut_
49  << " Set minimum noise value to: " << minStripNoiseCut_
50  << " Set maximum noise value to: " << maxStripNoiseCut_
51  << " Set maximum noise significance value to: " << maxStripNoiseSignificanceCut_
52  << " Set minimum Anderson-Darling p-value to: " << adProbabCut_
53  << " Set minimum Kolmogorov-Smirnov p-value to: " << ksProbabCut_
54  << " Set minimum Jacque-Bera p-value to: " << jbProbabCut_
55  << " Set minimum Chi2 p-value to: " << chi2ProbabCut_
56  << " Set N-sigma for the integral to : " << integralNsigma_
57  << " Set maximum integral tail at N-sigma to : " << integralTailCut_
58  << " Set maximum Kurtosis to : " << kurtosisCut_;
59 }
60 
61 // ----------------------------------------------------------------------------
62 //
63 
64 void PedsFullNoiseAlgorithm::extract( const std::vector<TH1*>& histos ) {
65 
66  if ( !anal() ) {
67  edm::LogWarning(mlCommissioning_)
68  << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
69  << " NULL pointer to Analysis object!";
70  return;
71  }
72 
73  // Check number of histograms --> Pedestal, noise and noise2D
74  if ( histos.size() != 3 ) {
75  anal()->addErrorCode(sistrip::numberOfHistos_);
76  }
77 
78  // Extract FED key from histo title --> i.e. APV pairs or LLD channel
79  if ( !histos.empty() ) {
80  anal()->fedKey( extractFedKey( histos.front() ) );
81  }
82 
83  // Extract 1D histograms
84  std::vector<TH1*>::const_iterator ihis = histos.begin();
85  for ( ; ihis != histos.end(); ihis++ ) {
86 
87  // Check for NULL pointer
88  if ( !(*ihis) ) { continue; }
89 
90  SiStripHistoTitle title( (*ihis)->GetName() );
91  if ( title.runType() != sistrip::PEDS_FULL_NOISE ) {
92  anal()->addErrorCode(sistrip::unexpectedTask_);
93  continue;
94  }
95 
96  // Extract peds histos
97  if ( title.extraInfo().find(sistrip::extrainfo::roughPedestals_) != std::string::npos ) {
98  //@@ something here for rough peds?
99  }
100  else if ( title.extraInfo().find(sistrip::extrainfo::pedestals_) != std::string::npos ) {
101  hPeds_.first = *ihis;
102  hPeds_.second = (*ihis)->GetName();
103  }
104  else if ( title.extraInfo().find(sistrip::extrainfo::commonMode_) != std::string::npos ) {
105  //@@ something here for CM plots?
106  }
107  else if ( title.extraInfo().find(sistrip::extrainfo::noiseProfile_) != std::string::npos ) {
108  //@@ something here for noise profile plot?
109  hNoise_.first = *ihis;
110  hNoise_.second = (*ihis)->GetName();
111  }
112  else if ( title.extraInfo().find(sistrip::extrainfo::noise2D_) != std::string::npos ) {
113  hNoise2D_.first = *ihis;
114  hNoise2D_.second = (*ihis)->GetName();
115  }
116  else {
117  anal()->addErrorCode(sistrip::unexpectedExtraInfo_);
118  }
119  }
120 }
121 
122 // resetting vectors
123 void PedsFullNoiseAlgorithm::reset(PedsFullNoiseAnalysis* ana){
124 
125  for(size_t iapv = 0 ; iapv < ana->peds_.size(); iapv++){
126  ana->pedsMean_[iapv] = 0.;
127  ana->rawMean_[iapv] = 0.;
128  ana->noiseMean_[iapv] = 0.;
129  ana->pedsSpread_[iapv] = 0.;
130  ana->noiseSpread_[iapv] = 0.;
131  ana->rawSpread_[iapv] = 0.;
132  ana->pedsMax_[iapv] = 0.;
133  ana->pedsMin_[iapv] = 0.;
134  ana->rawMax_[iapv] = 0.;
135  ana->rawMin_[iapv] = 0.;
136  ana->noiseMax_[iapv] = 0.;
137  ana->noiseMin_[iapv] = 0.;
138 
139  for(size_t istrip = 0; istrip < ana->peds_[iapv].size(); istrip++){
140  ana->peds_[iapv][istrip] = 0.;
141  ana->noise_[iapv][istrip] = 0.;
142  ana->raw_[iapv][istrip] = 0.;
143  ana->adProbab_[iapv][istrip] = 0.;
144  ana->ksProbab_[iapv][istrip] = 0.;
145  ana->jbProbab_[iapv][istrip] = 0.;
146  ana->chi2Probab_[iapv][istrip] = 0.;
147  ana->residualRMS_[iapv][istrip] = 0.;
148  ana->residualSigmaGaus_[iapv][istrip] = 0.;
149  ana->noiseSignificance_[iapv][istrip] = 0.;
150  ana->residualMean_[iapv][istrip] = 0.;
151  ana->residualSkewness_[iapv][istrip] = 0.;
152  ana->residualKurtosis_[iapv][istrip] = 0.;
153  ana->residualIntegralNsigma_[iapv][istrip] = 0.;
154  ana->residualIntegral_[iapv][istrip] = 0.;
155  ana->deadStripBit_[iapv][istrip] = 0;
156  ana->badStripBit_[iapv][istrip] = 0;
157  }
158  }
159 }
160 
161 // -----------------------------------------------------------------------------
162 //
163 void PedsFullNoiseAlgorithm::analyse() {
164 
165  // check base analysis object
166  if ( !anal() ) {
167  edm::LogWarning(mlCommissioning_)
168  << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
169  << " NULL pointer to base Analysis object!";
170  return;
171  }
172 
173  CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>( anal() );
174  PedsFullNoiseAnalysis* ana = dynamic_cast<PedsFullNoiseAnalysis*>( tmp );
175 
176  // check PedsFullNoiseAnalysis object
177  if ( !ana ) {
178  edm::LogWarning(mlCommissioning_)
179  << "[PedsFullNoiseAlgorithm::" << __func__ << "]"
180  << " NULL pointer to derived Analysis object!";
181  return;
182  }
183 
184 
185  // check if the histograms exists
186  if ( !hPeds_.first) {
187  ana->addErrorCode(sistrip::nullPtr_);
188  return;
189  }
190 
191  if ( !hNoise_.first ) {
192  ana->addErrorCode(sistrip::nullPtr_);
193  return;
194  }
195 
196  if ( !hNoise2D_.first ) {
197  ana->addErrorCode(sistrip::nullPtr_);
198  return;
199  }
200 
201  // take the histograms
202  TProfile *histoPeds = dynamic_cast<TProfile *>(hPeds_.first);
203  TProfile *histoNoiseMean = dynamic_cast<TProfile *>(hNoise_.first);
204  TH2S * histoNoise = dynamic_cast<TH2S*>(hNoise2D_.first);
205 
206  // Make sanity checks about pointers
207  if (not histoPeds) {
208  ana->addErrorCode(sistrip::nullPtr_);
209  return;
210  }
211 
212  if (not histoNoiseMean) {
213  ana->addErrorCode(sistrip::nullPtr_);
214  return;
215  }
216 
217  if (not histoNoise) {
218  ana->addErrorCode(sistrip::nullPtr_);
219  return;
220  }
221 
222  // check the binning --> each x-axis bin is 1 strip -> 2APV per lldChannel -> 256 strips
223  if (histoPeds->GetNbinsX() != 256 ) {
224  ana->addErrorCode(sistrip::numberOfBins_);
225  return;
226  }
227 
228  //check the binning --> each x-axis bin is 1 strip -> 2APV per lldChannel -> 256 strips
229  if (histoNoiseMean->GetNbinsX() != 256 ) {
230  ana->addErrorCode(sistrip::numberOfBins_);
231  return;
232  }
233 
234  //check the binning --> each y-axis bin is 1 strip -> 2APV per lldChannel -> 256 strips
235  if (histoNoise->GetNbinsY() != 256 ) {
236  ana->addErrorCode(sistrip::numberOfBins_);
237  return;
238  }
239 
240  //Reset values
241  reset(ana);
242 
243 
244  // loop on each strip
245  uint32_t apvID = -1;
246 
247  // Save basic information at strip / APV level
248  vector<float> ped_max;
249  vector<float> ped_min;
250  vector<float> raw_max;
251  vector<float> raw_min;
252  vector<float> noise_max;
253  vector<float> noise_min;
254 
255  // loop on each strip in the lldChannel
256  for(int iStrip = 0; iStrip < histoPeds->GetNbinsX(); iStrip++){
257 
258  if(iStrip < histoPeds->GetNbinsX()/2) apvID = 0;
259  else apvID = 1;
260 
261  int stripBin = 0;
262  if(iStrip >= 128) stripBin = iStrip-128;
263  else stripBin = iStrip;
264 
265  ana->peds_[apvID][stripBin] = histoPeds->GetBinContent(iStrip+1); // pedestal value
266  ana->noise_[apvID][stripBin] = histoNoiseMean->GetBinContent(iStrip+1); // noise value
267  ana->raw_[apvID][stripBin] = histoPeds->GetBinError(iStrip+1); // raw noise value
268 
269  ana->pedsMean_[apvID] += ana->peds_[apvID][stripBin]; // mean pedestal
270  ana->rawMean_[apvID] += ana->raw_[apvID][stripBin]; // mean raw noise
271  ana->noiseMean_[apvID] += ana->noise_[apvID][stripBin]; // mean noise
272 
273  // max pedestal
274  if(ped_max.size() < apvID+1)
275  ped_max.push_back(ana->peds_[apvID][stripBin]);
276  else{
277  if(ana->peds_[apvID][stripBin] > ped_max.at(apvID))
278  ped_max.at(apvID) = ana->peds_[apvID][stripBin];
279  }
280 
281  // min pedestal
282  if(ped_min.size() < apvID+1)
283  ped_min.push_back(ana->peds_[apvID][stripBin]);
284  else{
285  if(ana->peds_[apvID][stripBin] < ped_min.at(apvID))
286  ped_min.at(apvID) = ana->peds_[apvID][stripBin]; // min pedestal
287  }
288 
289  // max noise
290  if(noise_max.size() < apvID+1)
291  noise_max.push_back(ana->noise_[apvID][stripBin]);
292  else{
293  if(ana->noise_[apvID][stripBin] > noise_max.at(apvID))
294  noise_max.at(apvID) = ana->noise_[apvID][stripBin];
295  }
296 
297  // min noise
298  if(noise_min.size() < apvID+1)
299  noise_min.push_back(ana->noise_[apvID][stripBin]);
300  else{
301  if(ana->noise_[apvID][stripBin] < noise_min.at(apvID))
302  noise_min.at(apvID) = ana->noise_[apvID][stripBin];
303  }
304 
305  // max raw
306  if(raw_max.size() < apvID+1)
307  raw_max.push_back(ana->raw_[apvID][stripBin]);
308  else{
309  if(ana->raw_[apvID][stripBin] > raw_max.at(apvID))
310  raw_max.at(apvID) = ana->raw_[apvID][stripBin];
311  }
312 
313  // min raw
314  if(raw_min.size() < apvID+1)
315  raw_min.push_back(ana->raw_[apvID][stripBin]);
316  else{
317  if(ana->raw_[apvID][stripBin] < raw_min.at(apvID))
318  raw_min.at(apvID) = ana->raw_[apvID][stripBin];
319  }
320  }
321 
322  // Mean values
323  for(unsigned int iApv = 0; iApv < ana->pedsMean_.size(); iApv++){
324  ana->pedsMean_.at(iApv) /= (ana->peds_[iApv].size()); // calculate mean pedestal per APV
325  ana->rawMean_.at(iApv) /= (ana->raw_[iApv].size()); // calculate mean raw noise per APV
326  ana->noiseMean_.at(iApv) /= (ana->noise_[iApv].size()); // calculate mean noise per APV
327  }
328 
329  // Min and Max
330  for(unsigned int iApv = 0; iApv < ped_max.size(); iApv++){
331  if(ped_max.at(iApv) > sistrip::maximum_)
332  ana->pedsMax_.at(iApv) = sistrip::maximum_;
333  else if(ped_max.at(iApv) < -1.*sistrip::maximum_)
334  ana->pedsMax_.at(iApv) = -1.*sistrip::maximum_;
335  else
336  ana->pedsMax_.at(iApv) = ped_max.at(iApv);
337 
338  if(ped_min.at(iApv) > sistrip::maximum_)
339  ana->pedsMin_.at(iApv) = sistrip::maximum_;
340  else if(ped_min.at(iApv) < -1.*sistrip::maximum_)
341  ana->pedsMin_.at(iApv) = -1.*sistrip::maximum_;
342  else
343  ana->pedsMin_.at(iApv) = ped_min.at(iApv);
344 
345  if(noise_max.at(iApv) > sistrip::maximum_)
346  ana->noiseMax_.at(iApv) = sistrip::maximum_;
347  else if(noise_max.at(iApv) < -1.*sistrip::maximum_)
348  ana->noiseMax_.at(iApv) = -1.*sistrip::maximum_;
349  else
350  ana->noiseMax_.at(iApv) = noise_max.at(iApv);
351 
352  if(noise_min.at(iApv) > sistrip::maximum_)
353  ana->noiseMin_.at(iApv) = sistrip::maximum_;
354  else if(noise_min.at(iApv) < -1.*sistrip::maximum_)
355  ana->noiseMin_.at(iApv) = -1.*sistrip::maximum_;
356  else
357  ana->noiseMin_.at(iApv) = noise_min.at(iApv);
358 
359  if(raw_max.at(iApv) > sistrip::maximum_)
360  ana->rawMax_.at(iApv) = sistrip::maximum_;
361  else if(raw_max.at(iApv) < -1.*sistrip::maximum_)
362  ana->rawMax_.at(iApv) = -1.*sistrip::maximum_;
363  else
364  ana->rawMax_.at(iApv) = raw_max.at(iApv);
365 
366  if(raw_min.at(iApv) > sistrip::maximum_)
367  ana->rawMin_.at(iApv) = sistrip::maximum_;
368  else if(raw_min.at(iApv) < -1.*sistrip::maximum_)
369  ana->rawMin_.at(iApv) = -1.*sistrip::maximum_;
370  else
371  ana->rawMin_.at(iApv) = raw_min.at(iApv);
372 
373  }
374 
375  // Calculate the spread for noise and pedestal
376  apvID = -1;
377 
378  for(int iStrip = 0; iStrip < histoNoiseMean->GetNbinsX(); iStrip++){
379  if(iStrip < histoNoiseMean->GetNbinsX()/2) apvID = 0;
380  else apvID = 1;
381  ana->pedsSpread_[apvID] += pow(histoPeds->GetBinContent(iStrip+1)-ana->pedsMean_.at(apvID),2);
382  ana->noiseSpread_[apvID] += pow(histoNoiseMean->GetBinContent(iStrip+1)-ana->noiseMean_.at(apvID),2);
383  ana->rawSpread_[apvID] += pow(histoPeds->GetBinError(iStrip+1)-ana->rawMean_.at(apvID),2);
384  }
385 
386  for(unsigned int iApv = 0; iApv < ana->pedsSpread_.size(); iApv++){
387  ana->pedsSpread_[iApv] = sqrt(ana->pedsSpread_[iApv])/sqrt(ana->peds_[iApv].size() -1);
388  ana->noiseSpread_[iApv] = sqrt(ana->noiseSpread_[iApv])/sqrt(ana->noise_[iApv].size()-1);
389  ana->rawSpread_[iApv] = sqrt(ana->rawSpread_[iApv])/sqrt(ana->raw_[iApv].size() -1);
390  }
391 
392  // loop on each strip in the lldChannel
393  apvID = 0;
394  TH1S* histoResidualStrip = new TH1S("histoResidualStrip","",histoNoise->GetNbinsX(),histoNoise->GetXaxis()->GetXmin(),histoNoise->GetXaxis()->GetXmax());
395  histoResidualStrip->Sumw2();
396  histoResidualStrip->SetDirectory(nullptr);
397  TF1* fitFunc = new TF1 ("fitFunc","gaus(0)",histoNoise->GetXaxis()->GetXmin(),histoNoise->GetXaxis()->GetXmax());
398  TF1* fit2Gaus = nullptr;
399  TH1F* randomHisto = nullptr;
400  TFitResultPtr result;
401 
402  for(int iStrip = 0; iStrip < histoNoise->GetNbinsY(); iStrip++){
403  // tell which APV
404  if(iStrip < histoNoise->GetNbinsY()/2) apvID = 0;
405  else apvID = 1;
406  histoResidualStrip->Reset();
407 
408  int stripBin = 0;
409  if(iStrip >= 128) stripBin = iStrip-128;
410  else stripBin = iStrip;
411 
412  for(int iBinX = 0; iBinX < histoNoise->GetNbinsX(); iBinX++){
413  histoResidualStrip->SetBinContent(iBinX+1,histoNoise->GetBinContent(iBinX+1,iStrip+1));
414  histoResidualStrip->SetBinError(iBinX+1,histoNoise->GetBinError(iBinX+1,iStrip+1));
415  }
416 
417  if(histoResidualStrip->Integral() == 0){ // dead strip --> i.e. no data
418 
419  // set default values
420  ana->adProbab_[apvID][stripBin] = 0;
421  ana->ksProbab_[apvID][stripBin] = 0;
422  ana->jbProbab_[apvID][stripBin] = 0;
423  ana->chi2Probab_[apvID][stripBin] = 0;
424  ana->noiseSignificance_[apvID][stripBin] = 0;
425  ana->residualMean_[apvID][stripBin] = 0;
426  ana->residualRMS_[apvID][stripBin] = 0;
427  ana->residualSigmaGaus_[apvID][stripBin] = 0;
428  ana->residualSkewness_[apvID][stripBin] = 0;
429  ana->residualKurtosis_[apvID][stripBin] = 0;
430  ana->residualIntegralNsigma_[apvID][stripBin] = 0;
431  ana->residualIntegral_[apvID][stripBin] = 0;
432  ana->deadStrip_[apvID].push_back(stripBin);
433  ana->deadStripBit_[apvID][stripBin] = 1;
434  ana->badStripBit_[apvID][stripBin] = 0;
435 
436  SiStripFecKey fec_key(ana->fecKey());
437  LogTrace(mlDqmClient_)<<"DeadStrip: fecCrate "
438  <<" "<<fec_key.fecCrate()
439  <<" fecSlot "<<fec_key.fecSlot()
440  <<" fecRing "<<fec_key.fecRing()
441  <<" ccuAddr "<<fec_key.ccuAddr()
442  <<" ccChan "<<fec_key.ccuChan()
443  <<" lldChan "<<fec_key.lldChan()
444  <<" apvID "<<apvID
445  <<" stripID "<<iStrip;
446 
447 
448  continue;
449  }
450 
451  // set / calculated basic quantities
452  ana->residualMean_[apvID][stripBin] = histoResidualStrip->GetMean();
453  ana->residualRMS_[apvID][stripBin] = histoResidualStrip->GetRMS();
454  ana->residualSkewness_[apvID][stripBin] = histoResidualStrip->GetSkewness();
455  ana->residualKurtosis_[apvID][stripBin] = histoResidualStrip->GetKurtosis();
456  ana->noiseSignificance_[apvID][stripBin] = (ana->noise_[apvID][stripBin]-ana->noiseMean_[apvID])/ana->noiseSpread_[apvID];
457  ana->residualIntegral_[apvID][stripBin] = histoResidualStrip->Integral();
458  ana->residualIntegralNsigma_[apvID][stripBin] =
459  (histoResidualStrip->Integral(histoResidualStrip->FindBin(ana->residualMean_[apvID][stripBin]+ana->residualRMS_[apvID][stripBin]*integralNsigma_),histoResidualStrip->GetNbinsX()+1) +
460  histoResidualStrip->Integral(0,histoResidualStrip->FindBin(ana->residualMean_[apvID][stripBin]-ana->residualRMS_[apvID][stripBin]*integralNsigma_)))/ana->residualIntegral_[apvID][stripBin];
461 
462  // performing a Gaussian fit of the residual distribution
463  fitFunc->SetRange(histoNoise->GetXaxis()->GetXmin(),histoNoise->GetXaxis()->GetXmax());
464  fitFunc->SetParameters(ana->residualIntegral_[apvID][stripBin],ana->residualMean_[apvID][stripBin],ana->residualRMS_[apvID][stripBin]);
465  result = histoResidualStrip->Fit(fitFunc,"QSRN");
466 
467  // Good gaussian fit
468  if(result.Get()){
469 
470  ana->residualSigmaGaus_[apvID][stripBin] = fitFunc->GetParameter(2);
471  ana->chi2Probab_[apvID][stripBin] = result->Prob();
472 
473  // jacque bera probability
474  float jbVal = (ana->residualIntegral_[apvID][stripBin]/6)*(pow(ana->residualSkewness_[apvID][stripBin],2)+pow(ana->residualKurtosis_[apvID][stripBin],2)/4);
475  ana->jbProbab_[apvID][stripBin] = ROOT::Math::chisquared_cdf_c(jbVal,2);
476 
477  //Kolmogorov Smirnov and Anderson Darlong
478  if(randomHisto == nullptr)
479  randomHisto = (TH1F*) histoResidualStrip->Clone("randomHisto");
480  randomHisto->Reset();
481  randomHisto->SetDirectory(nullptr);
482 
483  if(generateRandomHisto_){
484  randomHisto->FillRandom("fitFunc",histoResidualStrip->Integral());
485  if(randomHisto->Integral() != 0){
486  ana->ksProbab_[apvID][stripBin] = histoResidualStrip->KolmogorovTest(randomHisto,"N");
487  ana->adProbab_[apvID][stripBin] = histoResidualStrip->AndersonDarlingTest(randomHisto);
488  }
489  else{
490  ana->ksProbab_[apvID][stripBin] = 0;
491  ana->adProbab_[apvID][stripBin] = 0;
492  }
493 
494  }
495  else{
496  randomHisto->Add(fitFunc);
497  if(randomHisto->Integral() != 0){
498  ana->ksProbab_[apvID][stripBin] = histoResidualStrip->KolmogorovTest(randomHisto,"N");
499  ROOT::Fit::BinData data1;
500  ROOT::Fit::BinData data2;
501  ROOT::Fit::FillData(data1,histoResidualStrip,nullptr);
502  data2.Initialize(randomHisto->GetNbinsX()+1,1);
503  for(int ibin = 0; ibin < randomHisto->GetNbinsX(); ibin++){
504  if(histoResidualStrip->GetBinContent(ibin+1) != 0 or randomHisto->GetBinContent(ibin+1) >= 1)
505  data2.Add(randomHisto->GetBinCenter(ibin+1),randomHisto->GetBinContent(ibin+1),randomHisto->GetBinError(ibin+1));
506  }
507 
508  double probab, value;
509  ROOT::Math::GoFTest::AndersonDarling2SamplesTest(data1,data2,probab,value);
510  ana->adProbab_[apvID][stripBin] = probab;
511  }
512  else{
513  ana->ksProbab_[apvID][stripBin] = 0;
514  ana->adProbab_[apvID][stripBin] = 0;
515  }
516  }
517  }
518 
519  // start applying selections storing output
520  bool badStripFlag = false;
521  ana->deadStripBit_[apvID][stripBin] = 0; // is not dead if the strip has data
522 
523  if(fabs(ana->residualMean_[apvID][stripBin]) > maxDriftResidualCut_ and not badStripFlag) {//mean value
524  ana->shiftedStrip_[apvID].push_back(stripBin);
525  badStripFlag = true;
526  }
527 
528  if(ana->residualRMS_[apvID][stripBin] < minStripNoiseCut_ and not badStripFlag){ // low noise
529  ana->lowNoiseStrip_[apvID].push_back(stripBin);
530  badStripFlag = true;
531  }
532 
533  if(ana->residualRMS_[apvID][stripBin] > maxStripNoiseCut_ and not badStripFlag){ // large noise
534  ana->largeNoiseStrip_[apvID].push_back(stripBin);
535  badStripFlag = true;
536  }
537 
538  if(fabs(ana->noiseSignificance_[apvID][stripBin]) > maxStripNoiseSignificanceCut_ and not badStripFlag){ // large noise significance
539  ana->largeNoiseSignificance_[apvID].push_back(stripBin);
540  badStripFlag = true;
541  }
542 
543  if(result.Get() and result->Status() != 0) // bad fit status
544  ana->badFitStatus_[apvID].push_back(stripBin);
545 
546  if(ana->adProbab_[apvID][stripBin] < adProbabCut_ and not badStripFlag) // bad AD p-value --> store the strip-id
547  ana->badADProbab_[apvID].push_back(stripBin);
548 
549  if(ana->ksProbab_[apvID][stripBin] < ksProbabCut_ and not badStripFlag) // bad KS p-value --> store the strip-id
550  ana->badKSProbab_[apvID].push_back(stripBin);
551 
552  if(ana->jbProbab_[apvID][stripBin] < jbProbabCut_ and not badStripFlag) // bad JB p-value --> store the strip-id
553  ana->badJBProbab_[apvID].push_back(stripBin);
554 
555  if(ana->chi2Probab_[apvID][stripBin] < chi2ProbabCut_ and not badStripFlag) // bad CHI2 p-value --> store the strip-id
556  ana->badChi2Probab_[apvID].push_back(stripBin);
557 
558  if(ana->adProbab_[apvID][stripBin] < adProbabCut_ and ana->ksProbab_[apvID][stripBin] < ksProbabCut_ and
559  ana->jbProbab_[apvID][stripBin] < jbProbabCut_ and ana->chi2Probab_[apvID][stripBin] < chi2ProbabCut_)
560  badStripFlag = true; // bad strip is flagged as bad by all the methods
561 
562  if(ana->residualKurtosis_[apvID][stripBin] > kurtosisCut_ and ana->residualIntegralNsigma_[apvID][stripBin] > integralTailCut_ and not badStripFlag){ // bad tails
563  ana->badTailStrip_[apvID].push_back(stripBin);
564  badStripFlag = true;
565  }
566 
567  if(badStripFlag){ // loop for double peaked
568 
569  fit2Gaus = new TF1("dgaus","[0]*exp(-((x-[1])*(x-[1]))/(2*[2]*[2]))+[3]*exp(-((x-[4])*(x-[4]))/(2*[5]*[5]))",histoNoise->GetXaxis()->GetXmin(),histoNoise->GetXaxis()->GetXmax());
570  fit2Gaus->SetParameter(0,fitFunc->GetParameter(0)/2);
571  fit2Gaus->SetParameter(3,fitFunc->GetParameter(0)/2);
572  fit2Gaus->SetParameter(1,1.);
573  fit2Gaus->SetParameter(4,-1.);
574  fit2Gaus->SetParameter(2,fitFunc->GetParameter(2));
575  fit2Gaus->SetParameter(5,fitFunc->GetParameter(2));
576  fit2Gaus->SetParLimits(1,0.,histoNoise->GetXaxis()->GetXmax());
577  fit2Gaus->SetParLimits(4,histoNoise->GetXaxis()->GetXmin(),0);
578  result = histoResidualStrip->Fit(fit2Gaus,"QSR");
579 
580  // ashman distance
581  float ashman = TMath::Power(2,0.5)*abs(fit2Gaus->GetParameter(1)-fit2Gaus->GetParameter(4))/(sqrt(pow(fit2Gaus->GetParameter(2),2)+pow(fit2Gaus->GetParameter(5),2)));
582  // amplitude
583  float amplitudeRatio = std::min(fit2Gaus->GetParameter(0),fit2Gaus->GetParameter(3))/std::max(fit2Gaus->GetParameter(0),fit2Gaus->GetParameter(3));
584 
585  if(ashman > ashmanDistance_ and amplitudeRatio > amplitudeRatio_)
586  ana->badDoublePeakStrip_[apvID].push_back(stripBin);
587  }
588 
589  if(badStripFlag){ // save the final bit
590 
591  ana->badStrip_[apvID].push_back(stripBin);
592  ana->badStripBit_[apvID][stripBin] = 1;
593 
594  SiStripFecKey fec_key(ana->fecKey());
595  LogTrace(mlDqmClient_)<<"BadStrip: fecCrate "
596  <<" "<<fec_key.fecCrate()
597  <<" fecSlot "<<fec_key.fecSlot()
598  <<" fecRing "<<fec_key.fecRing()
599  <<" ccuAddr "<<fec_key.ccuAddr()
600  <<" ccChan "<<fec_key.ccuChan()
601  <<" lldChan "<<fec_key.lldChan()
602  <<" apvID "<<apvID
603  <<" stripID "<<stripBin;
604 
605  }
606  else
607  ana->badStripBit_[apvID][stripBin] = 0;
608  }
609 
610 
611  ped_max.clear();
612  ped_min.clear();
613  raw_max.clear();
614  raw_min.clear();
615  noise_max.clear();
616  noise_min.clear();
617  if(histoResidualStrip) delete histoResidualStrip;
618  if(fitFunc) delete fitFunc;
619  if(randomHisto) delete randomHisto;
620  if(fit2Gaus) delete fit2Gaus;
621 
622 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:670
PedsFullNoiseAnalysis::peds_
VVFloat peds_
Quantitles that are always filled for every strip.
Definition: PedsFullNoiseAnalysis.h:107
sistrip::unexpectedTask_
static const char unexpectedTask_[]
Definition: ConstantsForCommissioningAnalysis.h:21
CommissioningAnalysis::fedKey
const uint32_t & fedKey() const
Definition: CommissioningAnalysis.h:141
mps_fire.result
result
Definition: mps_fire.py:291
SiStripHistoTitle
Utility class that holds histogram title.
Definition: SiStripHistoTitle.h:20
sistrip::numberOfHistos_
static const char numberOfHistos_[]
Definition: ConstantsForCommissioningAnalysis.h:16
nullptr
#define nullptr
Definition: GCC11Compatibility.h:37
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
sistrip::mlDqmClient_
static const char mlDqmClient_[]
Definition: ConstantsForLogger.h:19
sistrip::unexpectedExtraInfo_
static const char unexpectedExtraInfo_[]
Definition: ConstantsForCommissioningAnalysis.h:22
PedsFullNoiseAnalysis::deadStrip_
VVInt deadStrip_
Quantities filled only for bad strips i.e. vectors of strip-id.
Definition: PedsFullNoiseAnalysis.h:127
sistrip::numberOfBins_
static const char numberOfBins_[]
Definition: ConstantsForCommissioningAnalysis.h:18
sistrip
sistrip classes
Definition: EnsembleCalibrationLA.cc:11
PedsFullNoiseAnalysis
Histogram-based analysis for pedestal run.
Definition: PedsFullNoiseAnalysis.h:16
SiStripFecKey
Utility class that identifies a position within the strip tracker control structure, down to the level of an APV25.
Definition: SiStripFecKey.h:45
sistrip::extrainfo::noise2D_
static const char noise2D_[]
Definition: ConstantsForDqm.h:51
sistrip::mlCommissioning_
static const char mlCommissioning_[]
Definition: ConstantsForLogger.h:15
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
CommissioningAlgorithm::extractFedKey
uint32_t extractFedKey(const TH1 *const )
Definition: CommissioningAlgorithm.cc:31
or
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
sistrip::maximum_
static const uint16_t maximum_
Definition: Constants.h:20
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
sistrip::extrainfo::roughPedestals_
static const char roughPedestals_[]
Definition: ConstantsForDqm.h:49
CommissioningAnalysis::addErrorCode
virtual void addErrorCode(const std::string &error)
Definition: CommissioningAnalysis.h:155
CommissioningAnalysis::fecKey
const uint32_t & fecKey() const
Definition: CommissioningAnalysis.h:140
min
T min(T a, T b)
Definition: MathUtil.h:58
sistrip::extrainfo::commonMode_
static const char commonMode_[]
Definition: ConstantsForDqm.h:48
sistrip::extrainfo::noiseProfile_
static const char noiseProfile_[]
Definition: ConstantsForDqm.h:50
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:671
sistrip::extrainfo::pedestals_
static const char pedestals_[]
Definition: ConstantsForDqm.h:45
electrons_cff.bool
bool
Definition: electrons_cff.py:397
PedsFullNoiseAlgorithm::extract
void extract(const std::vector< TH1 * > &) override
Definition: PedsFullNoiseAlgorithm.cc:64
PedsFullNoiseAlgorithm::reset
void reset(PedsFullNoiseAnalysis *)
Definition: PedsFullNoiseAlgorithm.cc:123
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
CommissioningAnalysis
Abstract base for derived classes that provide analysis of commissioning histograms.
Definition: CommissioningAnalysis.h:18
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
CommissioningAlgorithm::anal
CommissioningAnalysis *const anal() const
Definition: CommissioningAlgorithm.h:55
sistrip::nullPtr_
static const char nullPtr_[]
Definition: ConstantsForCommissioningAnalysis.h:17
Generated for CMSSW Reference Manual by doxygen 1.8.11