CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Attributes
DTTimeBoxFitter Class Reference

#include <DTTimeBoxFitter.h>

Public Member Functions

 DTTimeBoxFitter (const TString &debugFileName=TString(""))
 Constructor. More...
 
std::pair< double, double > fitTimeBox (TH1F *hTimeBox)
 Fit the rising edge of the time box returning mean value and sigma (first and second respectively) More...
 
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. More...
 
void getInteractiveFitSeeds (TH1F *hTBox, double &mean, double &sigma, double &tBoxMax, double &xFitMin, double &xFitMax)
 Ask the user to provide the seeds. More...
 
void setFitSigma (double sigma)
 
void setInteractiveFit (bool isInteractive)
 Switch to interactive fit. More...
 
void setRebinning (int reb)
 Set the rebin. More...
 
void setVerbosity (unsigned int lvl)
 Set the verbosity of the output: 0 = silent, 1 = info, 2 = debug. More...
 
virtual ~DTTimeBoxFitter ()
 Destructor. More...
 

Private Attributes

TFile * hDebugFile
 
bool interactiveFit
 
int rebin
 
double theSigma
 
unsigned int theVerbosityLevel
 

Detailed Description

Fit the rising edge of the time box with the integral of a gaussian returning the mean value and the sigma.

Date:
2008/11/05 20:25:25
Revision:
1.7
Author
G. Cerminara - INFN Torino

Definition at line 20 of file DTTimeBoxFitter.h.

Constructor & Destructor Documentation

DTTimeBoxFitter::DTTimeBoxFitter ( const TString &  debugFileName = TString(""))

Constructor.

Definition at line 23 of file DTTimeBoxFitter.cc.

References hDebugFile, interactiveFit, and rebin.

23  : hDebugFile(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 }
unsigned int theVerbosityLevel
DTTimeBoxFitter::~DTTimeBoxFitter ( )
virtual

Destructor.

Definition at line 34 of file DTTimeBoxFitter.cc.

References hDebugFile.

34  {
35  if(hDebugFile != 0) hDebugFile->Close();
36 }

Member Function Documentation

pair< double, double > DTTimeBoxFitter::fitTimeBox ( TH1F *  hTimeBox)

Fit the rising edge of the time box returning mean value and sigma (first and second respectively)

Compute the ttrig (in ns) from the Time Box.

Definition at line 42 of file DTTimeBoxFitter.cc.

References gather_cfg::cout, getFitSeeds(), getInteractiveFitSeeds(), interactiveFit, intGauss(), timingPdfMaker::mean, rebin, and theVerbosityLevel.

42  {
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 }
double intGauss(double *x, double *par)
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.
void getInteractiveFitSeeds(TH1F *hTBox, double &mean, double &sigma, double &tBoxMax, double &xFitMin, double &xFitMax)
Ask the user to provide the seeds.
unsigned int theVerbosityLevel
tuple cout
Definition: gather_cfg.py:121
void DTTimeBoxFitter::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.

Definition at line 136 of file DTTimeBoxFitter.cc.

References abs, filterCSVwithJSON::copy, gather_cfg::cout, delta, diffTreeTool::diff, hDebugFile, i, MergeJob_cfg::interval, j, GetRecoTauVFromDQM_MC_cff::next, errorMatrix2Lands_multiChannel::start, theSigma, theVerbosityLevel, and dtDQMClient_cfg::threshold.

Referenced by fitTimeBox().

137  {
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 }
dbl * delta
Definition: mlp_gen.cc:36
int i
Definition: DBlmapReader.cc:9
tuple interval
Definition: MergeJob_cfg.py:20
#define abs(x)
Definition: mlp_lapack.h:159
int j
Definition: DBlmapReader.cc:9
unsigned int theVerbosityLevel
tuple cout
Definition: gather_cfg.py:121
void DTTimeBoxFitter::getInteractiveFitSeeds ( TH1F *  hTBox,
double &  mean,
double &  sigma,
double &  tBoxMax,
double &  xFitMin,
double &  xFitMax 
)

Ask the user to provide the seeds.

Definition at line 109 of file DTTimeBoxFitter.cc.

References gather_cfg::cout, timingPdfMaker::mean, theSigma, and theVerbosityLevel.

Referenced by fitTimeBox().

110  {
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 }
unsigned int theVerbosityLevel
tuple cout
Definition: gather_cfg.py:121
void DTTimeBoxFitter::setFitSigma ( double  sigma)
inline

Definition at line 58 of file DTTimeBoxFitter.h.

References theSigma.

58  {
59  theSigma = sigma;
60  }
void DTTimeBoxFitter::setInteractiveFit ( bool  isInteractive)
inline

Switch to interactive fit.

Definition at line 49 of file DTTimeBoxFitter.h.

References interactiveFit.

49  {
50  interactiveFit = isInteractive;
51  }
void DTTimeBoxFitter::setRebinning ( int  reb)
inline

Set the rebin.

Definition at line 54 of file DTTimeBoxFitter.h.

References rebin.

54  {
55  rebin = reb;
56  }
void DTTimeBoxFitter::setVerbosity ( unsigned int  lvl)
inline

Set the verbosity of the output: 0 = silent, 1 = info, 2 = debug.

Definition at line 44 of file DTTimeBoxFitter.h.

References theVerbosityLevel.

44  {
45  theVerbosityLevel = lvl;
46  }
unsigned int theVerbosityLevel

Member Data Documentation

TFile* DTTimeBoxFitter::hDebugFile
private

Definition at line 66 of file DTTimeBoxFitter.h.

Referenced by DTTimeBoxFitter(), getFitSeeds(), and ~DTTimeBoxFitter().

bool DTTimeBoxFitter::interactiveFit
private

Definition at line 69 of file DTTimeBoxFitter.h.

Referenced by DTTimeBoxFitter(), fitTimeBox(), and setInteractiveFit().

int DTTimeBoxFitter::rebin
private

Definition at line 70 of file DTTimeBoxFitter.h.

Referenced by DTTimeBoxFitter(), fitTimeBox(), and setRebinning().

double DTTimeBoxFitter::theSigma
private

Definition at line 71 of file DTTimeBoxFitter.h.

Referenced by getFitSeeds(), getInteractiveFitSeeds(), and setFitSigma().

unsigned int DTTimeBoxFitter::theVerbosityLevel
private

Definition at line 68 of file DTTimeBoxFitter.h.

Referenced by fitTimeBox(), getFitSeeds(), getInteractiveFitSeeds(), and setVerbosity().