CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DTTimeBoxFitter.cc
Go to the documentation of this file.
1 
2 /*
3  * See header file for a description of this class.
4  *
5  * $Date: 2010/01/08 12:51:42 $
6  * $Revision: 1.12 $
7  * \author G. Cerminara - INFN Torino
8  */
9 
11 
12 #include <iostream>
13 #include <vector>
14 
15 #include "TFile.h"
16 #include "TH1F.h"
17 #include "TMath.h"
18 #include "TF1.h"
19 #include "TString.h"
20 
21 using namespace std;
22 
23 DTTimeBoxFitter::DTTimeBoxFitter(const TString& debugFileName) : hDebugFile(0),
24  theVerbosityLevel(0),
25  theSigma(10.) {
26  // Create a root file for debug output only if needed
27  if(debugFileName != "") hDebugFile = new TFile(debugFileName.Data(), "RECREATE");
28  interactiveFit = false;
29  rebin = 1;
30 }
31 
32 
33 
35  if(hDebugFile != 0) hDebugFile->Close();
36 }
37 
38 
39 
40 
42 pair<double, double> DTTimeBoxFitter::fitTimeBox(TH1F *hTimeBox) {
43 
44  hTimeBox->Rebin(rebin);
45 
46  // Check if the histo contains any entry
47  if(hTimeBox->GetEntries() == 0) {
48  cout << "[DTTimeBoxFitter]***Error: the time box contains no entry!" << endl;
49  return make_pair(-1, -1);
50  }
51 
52  if(hTimeBox->GetEntries() < 50000) {
53  hTimeBox->Rebin(2);
54  }
55 
56  // Get seeds for the fit
57  // The TimeBox range to be fitted (the rising edge)
58  double xFitMin=0; // Min value for the fit
59  double xFitMax=0; // Max value for the fit
60  double xValue=0; // The mean value of the gaussian
61  double xFitSigma=0; // The sigma of the gaussian
62  double tBoxMax=0; // The max of the time box, it is used as seed for gaussian integral
63 
64  //hTimeBox->Rebin(2); //FIXME: Temporary for low statistics
65 
66  TH1F *hTimeBoxForSeed = (TH1F*) hTimeBox->Clone(); //FIXME: test
67  if(!interactiveFit) {
68  getFitSeeds(hTimeBoxForSeed, xValue, xFitSigma, tBoxMax, xFitMin, xFitMax);
69  } else {
70  getInteractiveFitSeeds(hTimeBoxForSeed, xValue, xFitSigma, tBoxMax, xFitMin, xFitMax);
71  }
72 
73  // Define the fitting function and use fit seeds
74  TF1 *fIntGaus = new TF1("IntGauss", intGauss, xFitMin, xFitMax, 3);
75  fIntGaus->SetParName(0, "Constant");
76  fIntGaus->SetParameter(0, tBoxMax);
77  fIntGaus->SetParName(1, "Mean");
78  fIntGaus->SetParameter(1, xValue);
79  fIntGaus->SetParName(2, "Sigma");
80  fIntGaus->SetParameter(2, xFitSigma);
81  fIntGaus->SetLineColor(kRed);
82 
83 
84  // Fit the histo
85  string option = "Q";
86  if(theVerbosityLevel >= 2)
87  option = "";
88 
89  hTimeBox->Fit("IntGauss", option.c_str(), "",xFitMin, xFitMax);
90 
91  // Get fitted parameters
92  double mean = fIntGaus->GetParameter("Mean");
93  double sigma = fIntGaus->GetParameter("Sigma");
94  // double constant = fIntGaus->GetParameter("Constant");
95  double chiSquare = fIntGaus->GetChisquare()/fIntGaus->GetNDF();
96 
97  if(theVerbosityLevel >= 1) {
98  cout << " === Fit Results: " << endl;
99  cout << " Fitted mean = " << mean << endl;
100  cout << " Fitted sigma = " << sigma << endl;
101  cout << " Reduced Chi Square = " << chiSquare << endl;
102  }
103  return make_pair(mean, sigma);
104 }
105 
106 
107 
108 // Get the seeds for the fit as input from user!It is used if interactiveFit == true
109 void DTTimeBoxFitter::getInteractiveFitSeeds(TH1F *hTBox, double& mean, double& sigma, double& tBoxMax,
110  double& xFitMin, double& xFitMax) {
111  if(theVerbosityLevel >= 1)
112  cout << " === Insert seeds for the Time Box fit:" << endl;
113 
114  cout << "Inser the fit mean:" << endl;
115  cin >> mean;
116 
117  sigma = theSigma; //FIXME: estimate it!
118 
119  tBoxMax = hTBox->GetMaximum();
120 
121  // Define the fit range
122  xFitMin = mean-5.*sigma;
123  xFitMax = mean+5.*sigma;
124 
125  if(theVerbosityLevel >= 1) {
126  cout << " Time Box Rising edge: " << mean << endl;
127  cout << " = Seeds and range for fit:" << endl;
128  cout << " Seed mean = " << mean << endl;
129  cout << " Seed sigma = " << sigma << endl;
130  cout << " Fitting from = " << xFitMin << " to " << xFitMax << endl << endl;
131  }
132 }
133 
134 
135 // Automatically compute the seeds the range to be used for time box fit
136 void DTTimeBoxFitter::getFitSeeds(TH1F *hTBox, double& mean, double& sigma, double& tBoxMax,
137  double& xFitMin, double& xFitMax) {
138  if(theVerbosityLevel >= 1)
139  cout << " === Looking for fit seeds in Time Box:" << endl;
140 
141 
142  // The approximate width of the time box
143  static const int tBoxWidth = 400; //FIXE: tune it
144 
145  int nBins = hTBox->GetNbinsX();
146  const int xMin = (int)hTBox->GetXaxis()->GetXmin();
147  const int xMax = (int)hTBox->GetXaxis()->GetXmax();
148  const int nEntries = (int)hTBox->GetEntries();
149 
150  double binValue = (double)(xMax-xMin)/(double)nBins;
151 
152  // Compute a threshold for TimeBox discrimination
153  const double threshold = binValue*nEntries/(double)(tBoxWidth*3.);
154  if(theVerbosityLevel >= 2)
155  cout << " Threshold for logic time box is (# entries): " << threshold << endl;
156 
157  int nRebins = 0; // protection for infinite loop
158  while(threshold > hTBox->GetMaximum()/2. && nRebins < 5) {
159  cout << " Rebinning!" << endl;
160  hTBox->Rebin(2);
161  nBins = hTBox->GetNbinsX();
162  binValue = (double)(xMax-xMin)/(double)nBins;
163  nRebins++;
164  }
165 
166  if(hDebugFile != 0) hDebugFile->cd();
167  TString hLName = TString(hTBox->GetName())+"L";
168  TH1F hLTB(hLName.Data(), "Logic Time Box", nBins, xMin, xMax);
169  // Loop over all time box bins and discriminate them accordigly to the threshold
170  for(int i = 1; i <= nBins; i++) {
171  if(hTBox->GetBinContent(i) > threshold)
172  hLTB.SetBinContent(i, 1);
173  }
174  if(hDebugFile != 0) hLTB.Write();
175 
176  // Look for the time box in the "logic histo" and save beginning and lenght of each plateau
177  vector< pair<int, int> > startAndLenght;
178  if(theVerbosityLevel >= 2)
179  cout << " Look for rising and folling edges of logic time box: " << endl;
180  int start = -1;
181  int lenght = -1;
182  for(int j = 1; j < nBins;j++) {
183  int diff = (int)hLTB.GetBinContent(j+1)-(int)hLTB.GetBinContent(j);
184  if(diff == 1) { // This is a rising edge
185  start = j;
186  lenght = 1;
187  if(theVerbosityLevel >= 2) {
188  cout << " >>>----" << endl;
189  cout << " bin: " << j << " is a rising edge" << endl;
190  }
191  } else if(diff == -1) { // This is a falling edge
192  startAndLenght.push_back(make_pair(start, lenght));
193  if(theVerbosityLevel >= 2) {
194  cout << " bin: " << j << " is a falling edge, lenght is: " << lenght << endl;
195  cout << " <<<----" << endl;
196  }
197  start = -1;
198  lenght = -1;
199  } else if(diff == 0 && (int)hLTB.GetBinContent(j) == 1) { // This is a bin within the plateau
200  lenght ++;
201  }
202  }
203 
204  if(theVerbosityLevel >= 2)
205  cout << " Add macro intervals: " << endl;
206  // Look for consecutive plateau: 2 plateau are consecutive if the gap is < 5 ns
207  vector< pair<int, int> > superIntervals;
208  for(vector< pair<int, int> >::const_iterator interval = startAndLenght.begin();
209  interval != startAndLenght.end();
210  ++interval) {
211  pair<int, int> theInterval = (*interval);
212  vector< pair<int, int> >::const_iterator next = interval;
213  while(++next != startAndLenght.end()) {
214  int gap = (*next).first - (theInterval.first+theInterval.second);
215  double gabInNs = binValue*gap;
216  if(theVerbosityLevel >= 2)
217  cout << " gap: " << gabInNs << "(ns)" << endl;
218  if(gabInNs > 20) {
219  break;
220  } else {
221  theInterval = make_pair(theInterval.first, theInterval.second+gap+(*next).second);
222  superIntervals.push_back(theInterval);
223  if(theVerbosityLevel >= 2)
224  cout << " Add interval, start: " << theInterval.first
225  << " width: " << theInterval.second << endl;
226  }
227  }
228  }
229 
230  // merge the vectors of intervals
231  copy(superIntervals.begin(), superIntervals.end(), back_inserter(startAndLenght));
232 
233  // Look for the plateau of the right lenght
234  if(theVerbosityLevel >= 2)
235  cout << " Look for the best interval:" << endl;
236  int delta = 999999;
237  int beginning = -1;
238  int tbWidth = -1;
239  for(vector< pair<int, int> >::const_iterator stAndL = startAndLenght.begin();
240  stAndL != startAndLenght.end();
241  stAndL++) {
242  if(abs((*stAndL).second - tBoxWidth) < delta) {
243  delta = abs((*stAndL).second - tBoxWidth);
244  beginning = (*stAndL).first;
245  tbWidth = (*stAndL).second;
246  if(theVerbosityLevel >= 2)
247  cout << " Candidate: First: " << beginning
248  << ", width: " << tbWidth
249  << ", delta: " << delta << endl;
250  }
251  }
252 
253  mean = xMin + beginning*binValue;
254  sigma = theSigma; //FIXME: estimate it!
255 
256  tBoxMax = hTBox->GetMaximum();
257 
258  // Define the fit range
259  xFitMin = mean-5.*sigma;
260  xFitMax = mean+5.*sigma;
261 
262  if(theVerbosityLevel >= 1) {
263  cout << " Time Box Rising edge: " << mean << endl;
264  cout << " Time Box Width: " << tbWidth*binValue << endl;
265  cout << " = Seeds and range for fit:" << endl;
266  cout << " Seed mean = " << mean << endl;
267  cout << " Seed sigma = " << sigma << endl;
268  cout << " Fitting from = " << xFitMin << " to " << xFitMax << endl << endl;
269  }
270 }
271 
272 
273 
274 double intGauss(double *x, double *par) {
275  double media = par[1];
276  double sigma = par[2];
277  double var = (x[0]-media)/(sigma*sqrt(2.));
278 
279  return 0.5*par[0]*(1+TMath::Erf(var));
280 
281 }
dbl * delta
Definition: mlp_gen.cc:36
int i
Definition: DBlmapReader.cc:9
tuple start
Check for commandline option errors.
Definition: dqm_diff.py:58
tuple interval
Definition: MergeJob_cfg.py:20
double intGauss(double *x, double *par)
#define abs(x)
Definition: mlp_lapack.h:159
T sqrt(T t)
Definition: SSEVec.h:48
int j
Definition: DBlmapReader.cc:9
void getFitSeeds(TH1F *hTBox, double &mean, double &sigma, double &tBoxMax, double &xFitMin, double &xFitMax)
Automatically compute the seeds the range to be used for time box fit.
virtual ~DTTimeBoxFitter()
Destructor.
void getInteractiveFitSeeds(TH1F *hTBox, double &mean, double &sigma, double &tBoxMax, double &xFitMin, double &xFitMax)
Ask the user to provide the seeds.
DTTimeBoxFitter(const TString &debugFileName=TString(""))
Constructor.
unsigned int theVerbosityLevel
std::pair< double, double > fitTimeBox(TH1F *hTimeBox)
Fit the rising edge of the time box returning mean value and sigma (first and second respectively) ...
tuple cout
Definition: gather_cfg.py:121
Definition: DDAxes.h:10