CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CalibrationAlgorithm.cc
Go to the documentation of this file.
7 #include "TProfile.h"
8 #include "TF1.h"
9 #include "TH1.h"
10 #include "TVirtualFitter.h"
11 #include <iostream>
12 #include <sstream>
13 #include <iomanip>
14 #include <cmath>
15 
16 using namespace sistrip;
17 
18 // ----------------------------------------------------------------------------
19 //
21  : CommissioningAlgorithm(anal),
22  deconv_fitter_(0),
23  peak_fitter_(0),
24  cal_(0)
25 {
26  deconv_fitter_ = new TF1("deconv_fitter",fdeconv_convoluted,-50,50,5);
27  deconv_fitter_->FixParameter(0,0);
28  deconv_fitter_->SetParLimits(1,-100,0);
29  deconv_fitter_->SetParLimits(2,0,200);
30  deconv_fitter_->SetParLimits(3,5,100);
31  deconv_fitter_->FixParameter(3,50);
32  deconv_fitter_->SetParLimits(4,0,50);
33  deconv_fitter_->SetParameters(0.,-10,0.96,50,20);
34  peak_fitter_ = new TF1("peak_fitter",fpeak_convoluted,-50,50,5);
35  peak_fitter_->FixParameter(0,0);
36  peak_fitter_->SetParLimits(1,-100,0);
37  peak_fitter_->SetParLimits(2,0,400);
38  peak_fitter_->SetParLimits(3,5,100);
39  peak_fitter_->FixParameter(3,50);
40  peak_fitter_->SetParLimits(4,0,50);
41  peak_fitter_->SetParameters(0.,-10,0.96,50,20);
42 }
43 
44 // ----------------------------------------------------------------------------
45 //
46 void CalibrationAlgorithm::extract( const std::vector<TH1*>& histos) {
47  if ( !anal() ) {
49  << "[CalibrationAlgorithm::" << __func__ << "]"
50  << " NULL pointer to base Analysis object!";
51  return;
52  }
53 
55  cal_ = dynamic_cast<CalibrationAnalysis*>( tmp );
56  if ( !cal_ ) {
58  << "[CalibrationAlgorithm::" << __func__ << "]"
59  << " NULL pointer to derived Analysis object!";
60  return;
61  }
62 
63  // Check number of histograms
64  if ( histos.size() != 32 && histos.size() !=2 ) {
66  }
67 
68  // Extract FED key from histo title
69  if ( !histos.empty() ) { cal_->fedKey( extractFedKey( histos.front() ) ); }
70 
71  // Extract histograms
72  std::vector<TH1*>::const_iterator ihis = histos.begin();
73  unsigned int cnt = 0;
74  for ( ; ihis != histos.end(); ihis++,cnt++ ) {
75 
76  // Check for NULL pointer
77  if ( !(*ihis) ) { continue; }
78 
79  // Check name
80  SiStripHistoTitle title( (*ihis)->GetName() );
81  if ( title.runType() != sistrip::CALIBRATION &&
82  title.runType() != sistrip::CALIBRATION_DECO &&
83  title.runType() != sistrip::CALIBRATION_SCAN &&
84  title.runType() != sistrip::CALIBRATION_SCAN_DECO ) {
86  continue;
87  }
88 
89  cal_->isScan_ = ( title.runType() == sistrip::CALIBRATION_SCAN ||
90  title.runType() == sistrip::CALIBRATION_SCAN_DECO );
91 
92  // Extract calibration histo
93  histo_[cnt].first = *ihis;
94  histo_[cnt].second = (*ihis)->GetTitle();
95  }
96 
97 }
98 
99 // ----------------------------------------------------------------------------
100 //
102 
103  if ( !cal_ ) {
105  << "[CalibrationAlgorithm::" << __func__ << "]"
106  << " NULL pointer to derived Analysis object!";
107  return;
108  }
109 
110  float Amean[2] = {0.,0.};
111  float Amin[2] = {2000.,2000.};
112  float Amax[2] = {0.,0.};
113  float Aspread[2] = {0.,0.};
114  float Tmean[2] = {0.,0.};
115  float Tmin[2] = {2000.,2000.};
116  float Tmax[2] = {0.,0.};
117  float Tspread[2] = {0.,0.};
118  float Rmean[2] = {0.,0.};
119  float Rmin[2] = {2000.,2000.};
120  float Rmax[2] = {0.,0.};
121  float Rspread[2] = {0.,0.};
122  float Cmean[2] = {0.,0.};
123  float Cmin[2] = {2000.,2000.};
124  float Cmax[2] = {0.,0.};
125  float Cspread[2] = {0.,0.};
126  float Smean[2] = {0.,0.};
127  float Smin[2] = {2000.,2000.};
128  float Smax[2] = {0.,0.};
129  float Sspread[2] = {0.,0.};
130  float Kmean[2] = {0.,0.};
131  float Kmin[2] = {2000000.,2000000.};
132  float Kmax[2] = {0.,0.};
133  float Kspread[2] = {0.,0.};
134 
135  unsigned int upperLimit = cal_->isScan_ ? 2 : 32;
136  float nStrips = cal_->isScan_ ? 1. : 16.;
137  for(unsigned int i=0;i<upperLimit;++i) {
138  if ( !histo_[i].first ) {
140  << " NULL pointer to histogram " << i << "!";
141  return;
142  }
143 
144  // determine which APV and strip is being looked at.
145  std::string title = histo_[i].first->GetName();
146  int apv = 0;
147  int strip = 0;
148  if(title.find("STRIP")!=std::string::npos && title.find("Apv")!=std::string::npos) {
149  strip = atoi(title.c_str()+title.find("STRIP")+6);
150  apv = (atoi(title.c_str()+title.find("Apv")+3))%2;
151  } else {
152  strip = (atoi(title.c_str()+title.find_last_of("_")+1))%16;
153  apv = (atoi(title.c_str()+title.find_last_of("_")+1))/16;
154  if(title.find("Apv")!=std::string::npos) {
155  apv = (atoi(title.c_str()+title.find("Apv")+3))%2;
156  strip = strip*8 + cal_->calchan_;
157  } else {
159  << " Malformed histogram title! Strip/APV not retreived: "
160  << title;
161  }
162  }
163 
164  // rescale the plot
165  correctDistribution(histo_[i].first);
166 
167  // amplitude
168  cal_->amplitude_[apv][strip] = histo_[i].first->GetMaximum();
169 
170  // rise time
171  cal_->riseTime_[apv][strip] = maximum(histo_[i].first) - turnOn(histo_[i].first);
172 
173  // tail 125 ns after the maximum
174  int lastBin = histo_[i].first->FindBin(histo_[i].first->GetBinCenter(histo_[i].first->GetMaximumBin())+125);
175  if(lastBin>histo_[i].first->GetNbinsX()-4) lastBin = histo_[i].first->GetNbinsX()-4;
176  if(histo_[i].first->GetMaximum()!=0)
177  cal_->tail_[apv][strip] = 100*histo_[i].first->GetBinContent(lastBin)/histo_[i].first->GetMaximum();
178  else
179  cal_->tail_[apv][strip] = 100;
180 
181  // perform the fit for the next quantities
182  TF1* fit = fitPulse(histo_[i].first);
183 
184  // time constant
185  cal_->timeConstant_[apv][strip] = fit->GetParameter(3);
186 
187  // smearing
188  cal_->smearing_[apv][strip] = fit->GetParameter(4);
189 
190  // chi2
191  cal_->chi2_[apv][strip] = fit->GetChisquare();
192 
193  //compute mean, max, min, spread
194  Amean[apv] += cal_->amplitude_[apv][strip]/nStrips;
195  Amin[apv] = Amin[apv]<cal_->amplitude_[apv][strip] ? Amin[apv] : cal_->amplitude_[apv][strip];
196  Amax[apv] = Amax[apv]>cal_->amplitude_[apv][strip] ? Amax[apv] : cal_->amplitude_[apv][strip];
197  Aspread[apv] += cal_->amplitude_[apv][strip]*cal_->amplitude_[apv][strip]/nStrips;
198  Tmean[apv] += cal_->tail_[apv][strip]/nStrips;
199  Tmin[apv] = Tmin[apv]<cal_->tail_[apv][strip] ? Tmin[apv] : cal_->tail_[apv][strip];
200  Tmax[apv] = Tmax[apv]>cal_->tail_[apv][strip] ? Tmax[apv] : cal_->tail_[apv][strip];
201  Tspread[apv] += cal_->tail_[apv][strip]*cal_->tail_[apv][strip]/nStrips;
202  Rmean[apv] += cal_->riseTime_[apv][strip]/nStrips;
203  Rmin[apv] = Rmin[apv]<cal_->riseTime_[apv][strip] ? Rmin[apv] : cal_->riseTime_[apv][strip];
204  Rmax[apv] = Rmax[apv]>cal_->riseTime_[apv][strip] ? Rmax[apv] : cal_->riseTime_[apv][strip];
205  Rspread[apv] += cal_->riseTime_[apv][strip]*cal_->riseTime_[apv][strip]/nStrips;
206  Cmean[apv] += cal_->timeConstant_[apv][strip]/nStrips;
207  Cmin[apv] = Cmin[apv]<cal_->timeConstant_[apv][strip] ? Cmin[apv] : cal_->timeConstant_[apv][strip];
208  Cmax[apv] = Cmax[apv]>cal_->timeConstant_[apv][strip] ? Cmax[apv] : cal_->timeConstant_[apv][strip];
209  Cspread[apv] += cal_->timeConstant_[apv][strip]*cal_->timeConstant_[apv][strip]/nStrips;
210  Smean[apv] += cal_->smearing_[apv][strip]/nStrips;
211  Smin[apv] = Smin[apv]<cal_->smearing_[apv][strip] ? Smin[apv] : cal_->smearing_[apv][strip];
212  Smax[apv] = Smax[apv]>cal_->smearing_[apv][strip] ? Smax[apv] : cal_->smearing_[apv][strip];
213  Sspread[apv] += cal_->smearing_[apv][strip]*cal_->smearing_[apv][strip]/nStrips;
214  Kmean[apv] += cal_->chi2_[apv][strip]/nStrips;
215  Kmin[apv] = Kmin[apv]<cal_->chi2_[apv][strip] ? Kmin[apv] : cal_->chi2_[apv][strip];
216  Kmax[apv] = Kmax[apv]>cal_->chi2_[apv][strip] ? Kmax[apv] : cal_->chi2_[apv][strip];
217  Kspread[apv] += cal_->chi2_[apv][strip]*cal_->chi2_[apv][strip]/nStrips;
218  }
219 
220  // fill the mean, max, min, spread, ... histograms.
221  for(int i=0;i<2;++i) {
222  cal_->mean_amplitude_[i] = Amean[i];
223  cal_->mean_tail_[i] = Tmean[i];
224  cal_->mean_riseTime_[i] = Rmean[i];
225  cal_->mean_timeConstant_[i] = Cmean[i];
226  cal_->mean_smearing_[i] = Smean[i];
227  cal_->mean_chi2_[i] = Kmean[i];
228  cal_->min_amplitude_[i] = Amin[i];
229  cal_->min_tail_[i] = Tmin[i];
230  cal_->min_riseTime_[i] = Rmin[i];
231  cal_->min_timeConstant_[i] = Cmin[i];
232  cal_->min_smearing_[i] = Smin[i];
233  cal_->min_chi2_[i] = Kmin[i];
234  cal_->max_amplitude_[i] = Amax[i];
235  cal_->max_tail_[i] = Tmax[i];
236  cal_->max_riseTime_[i] = Rmax[i];
237  cal_->max_timeConstant_[i] = Cmax[i];
238  cal_->max_smearing_[i] = Smax[i];
239  cal_->max_chi2_[i] = Kmax[i];
240  cal_->spread_amplitude_[i] = sqrt(fabs(Aspread[i]-Amean[i]*Amean[i]));
241  cal_->spread_tail_[i] = sqrt(fabs(Tspread[i]-Tmean[i]*Tmean[i]));
242  cal_->spread_riseTime_[i] = sqrt(fabs(Rspread[i]-Rmean[i]*Rmean[i]));
243  cal_->spread_timeConstant_[i] = sqrt(fabs(Cspread[i]-Cmean[i]*Cmean[i]));
244  cal_->spread_smearing_[i] = sqrt(fabs(Sspread[i]-Smean[i]*Smean[i]));
245  cal_->spread_chi2_[i] = sqrt(fabs(Kspread[i]-Kmean[i]*Kmean[i]));
246  }
247 }
248 
249 // ----------------------------------------------------------------------------
250 //
252 {
253  // return the curve
254  histo->Scale(-1);
255  if ( cal_ ) { if( cal_->isScan_ ) histo->Scale(1/16.); }
256 }
257 
258 // ----------------------------------------------------------------------------
259 //
261  float rangeLow,
262  float rangeHigh )
263 {
264  if(!cal_) return 0;
265  if(!histo) return 0;
266  float noise = 4.;
267  float N = round(histo->GetMaximum()/125.);
268  float error = sqrt(2*N)*noise;
269  //float error = sqrt(2)*noise*N; // ?????????
270  for(int i=1;i<=histo->GetNbinsX();++i) {
271  histo->SetBinError(i,error);
272  }
273  if (cal_->deconv_) {
274  if(rangeLow>rangeHigh)
275  histo->Fit(deconv_fitter_,"Q");
276  else
277  histo->Fit(deconv_fitter_,"0Q","",rangeLow,rangeHigh);
278  return deconv_fitter_;
279  } else {
280  if(rangeLow>rangeHigh)
281  histo->Fit(peak_fitter_,"Q");
282  else
283  histo->Fit(peak_fitter_,"0Q","",rangeLow,rangeHigh);
284  return peak_fitter_;
285  }
286 }
287 
288 // ----------------------------------------------------------------------------
289 //
291  int bin = h->GetMaximumBin();
292  // fit around the maximum with the detector response and take the max from the fit
293  TF1* fit = fitPulse(h,h->GetBinCenter(bin)-25,h->GetBinCenter(bin)+25);
294  return fit->GetMaximumX();
295 }
296 
297 // ----------------------------------------------------------------------------
298 //
300  // localize the rising edge
301  int bin=1;
302  float amplitude = h->GetMaximum();
303  for(; bin<= h->GetNbinsX() && h->GetBinContent(bin)<0.4*amplitude; ++bin) {}
304  float end = h->GetBinLowEdge(bin);
305  // fit the rising edge with a sigmoid
306  TF1* sigmoid = new TF1("sigmoid","[0]/(1+exp(-[1]*(x-[2])))+[3]",0,end);
307  sigmoid->SetParLimits(0,amplitude/10.,amplitude);
308  sigmoid->SetParLimits(1,0.05,0.5);
309  sigmoid->SetParLimits(2,end-10,end+10);
310  sigmoid->SetParLimits(3,-amplitude/10.,amplitude/10.);
311  sigmoid->SetParameters(amplitude/2.,0.1,end,0.);
312  h->Fit(sigmoid,"0QR");
313  // return the point where the fit = 3% signal.
314  float time = 0.;
315  float base = sigmoid->GetMinimum(0,end);
316  for(;time<end && (sigmoid->Eval(time)-base)<0.03*(amplitude-base);time += 0.1) {}
317  delete sigmoid;
318  return time-0.05;
319 }
320 
TF1 * fitPulse(TH1 *, float rangeLow=0, float rangeHigh=-1)
tuple base
Main Program
Definition: newFWLiteAna.py:92
int i
Definition: DBlmapReader.cc:9
static const char unexpectedTask_[]
double fdeconv_convoluted(double *x, double *par)
const uint32_t & fedKey() const
Utility class that holds histogram title.
static const char numberOfHistos_[]
CalibrationAnalysis * cal_
double fpeak_convoluted(double *x, double *par)
Analysis for calibration runs.
static const char mlCommissioning_[]
T sqrt(T t)
Definition: SSEVec.h:48
uint32_t extractFedKey(const TH1 *const )
virtual void addErrorCode(const std::string &error)
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
#define end
Definition: vmac.h:37
bool first
Definition: L1TdeRCT.cc:75
#define N
Definition: blowfish.cc:9
void correctDistribution(TH1 *) const
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
const double Rmax[kNumberCalorimeter]
virtual void extract(const std::vector< TH1 * > &)
Abstract base for derived classes that provide analysis of commissioning histograms.
const double Rmin[kNumberCalorimeter]
CommissioningAnalysis *const anal() const