CMS 3D CMS Logo

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.

Author
G. Cerminara - INFN Torino

Definition at line 18 of file DTTimeBoxFitter.h.

Constructor & Destructor Documentation

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

Constructor.

Definition at line 21 of file DTTimeBoxFitter.cc.

References hDebugFile, interactiveFit, and rebin.

21  : hDebugFile(0),
23  theSigma(10.) {
24  // Create a root file for debug output only if needed
25  if(debugFileName != "") hDebugFile = new TFile(debugFileName.Data(), "RECREATE");
26  interactiveFit = false;
27  rebin = 1;
28 }
unsigned int theVerbosityLevel
DTTimeBoxFitter::~DTTimeBoxFitter ( )
virtual

Destructor.

Definition at line 32 of file DTTimeBoxFitter.cc.

References hDebugFile.

32  {
33  if(hDebugFile != 0) hDebugFile->Close();
34 }

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 40 of file DTTimeBoxFitter.cc.

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

40  {
41 
42  hTimeBox->Rebin(rebin);
43 
44  // Check if the histo contains any entry
45  if(hTimeBox->GetEntries() == 0) {
46  cout << "[DTTimeBoxFitter]***Error: the time box contains no entry!" << endl;
47  return make_pair(-1, -1);
48  }
49 
50  if(hTimeBox->GetEntries() < 50000) {
51  hTimeBox->Rebin(2);
52  }
53 
54  // Get seeds for the fit
55  // The TimeBox range to be fitted (the rising edge)
56  double xFitMin=0; // Min value for the fit
57  double xFitMax=0; // Max value for the fit
58  double xValue=0; // The mean value of the gaussian
59  double xFitSigma=0; // The sigma of the gaussian
60  double tBoxMax=0; // The max of the time box, it is used as seed for gaussian integral
61 
62  //hTimeBox->Rebin(2); //FIXME: Temporary for low statistics
63 
64  TH1F *hTimeBoxForSeed = (TH1F*) hTimeBox->Clone(); //FIXME: test
65  if(!interactiveFit) {
66  getFitSeeds(hTimeBoxForSeed, xValue, xFitSigma, tBoxMax, xFitMin, xFitMax);
67  } else {
68  getInteractiveFitSeeds(hTimeBoxForSeed, xValue, xFitSigma, tBoxMax, xFitMin, xFitMax);
69  }
70 
71  // Define the fitting function and use fit seeds
72  TF1 *fIntGaus = new TF1("IntGauss", intGauss, xFitMin, xFitMax, 3);
73  fIntGaus->SetParName(0, "Constant");
74  fIntGaus->SetParameter(0, tBoxMax);
75  fIntGaus->SetParName(1, "Mean");
76  fIntGaus->SetParameter(1, xValue);
77  fIntGaus->SetParName(2, "Sigma");
78  fIntGaus->SetParameter(2, xFitSigma);
79  fIntGaus->SetLineColor(kRed);
80 
81 
82  // Fit the histo
83  string option = "Q";
84  if(theVerbosityLevel >= 2)
85  option = "";
86 
87  hTimeBox->Fit(fIntGaus, option.c_str(), "",xFitMin, xFitMax);
88 
89  // Get fitted parameters
90  double mean = fIntGaus->GetParameter("Mean");
91  double sigma = fIntGaus->GetParameter("Sigma");
92  // double constant = fIntGaus->GetParameter("Constant");
93  double chiSquare = fIntGaus->GetChisquare()/fIntGaus->GetNDF();
94 
95  if(theVerbosityLevel >= 1) {
96  cout << " === Fit Results: " << endl;
97  cout << " Fitted mean = " << mean << endl;
98  cout << " Fitted sigma = " << sigma << endl;
99  cout << " Reduced Chi Square = " << chiSquare << endl;
100  }
101  return make_pair(mean, sigma);
102 }
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
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 134 of file DTTimeBoxFitter.cc.

References funct::abs(), popcon2dropbox::copy(), gather_cfg::cout, delta, mps_update::diff, hDebugFile, mps_fire::i, createfilelist::int, GetRecoTauVFromDQM_MC_cff::next, theSigma, theVerbosityLevel, electronIdCutBased_cfi::threshold, anotherprimaryvertexanalyzer_cfi::xMax, and anotherprimaryvertexanalyzer_cfi::xMin.

Referenced by fitTimeBox().

135  {
136  if(theVerbosityLevel >= 1)
137  cout << " === Looking for fit seeds in Time Box:" << endl;
138 
139 
140  // The approximate width of the time box
141  static const int tBoxWidth = 400; //FIXE: tune it
142 
143  int nBins = hTBox->GetNbinsX();
144  const int xMin = (int)hTBox->GetXaxis()->GetXmin();
145  const int xMax = (int)hTBox->GetXaxis()->GetXmax();
146  const int nEntries = (int)hTBox->GetEntries();
147 
148  double binValue = (double)(xMax-xMin)/(double)nBins;
149 
150  // Compute a threshold for TimeBox discrimination
151  const double threshold = binValue*nEntries/(double)(tBoxWidth*3.);
152  if(theVerbosityLevel >= 2)
153  cout << " Threshold for logic time box is (# entries): " << threshold << endl;
154 
155  int nRebins = 0; // protection for infinite loop
156  while(threshold > hTBox->GetMaximum()/2. && nRebins < 5) {
157  cout << " Rebinning!" << endl;
158  hTBox->Rebin(2);
159  nBins = hTBox->GetNbinsX();
160  binValue = (double)(xMax-xMin)/(double)nBins;
161  nRebins++;
162  }
163 
164  if(hDebugFile != 0) hDebugFile->cd();
165  TString hLName = TString(hTBox->GetName())+"L";
166  TH1F hLTB(hLName.Data(), "Logic Time Box", nBins, xMin, xMax);
167  // Loop over all time box bins and discriminate them accordigly to the threshold
168  for(int i = 1; i <= nBins; i++) {
169  if(hTBox->GetBinContent(i) > threshold)
170  hLTB.SetBinContent(i, 1);
171  }
172  if(hDebugFile != 0) hLTB.Write();
173 
174  // Look for the time box in the "logic histo" and save beginning and lenght of each plateau
175  vector< pair<int, int> > startAndLenght;
176  if(theVerbosityLevel >= 2)
177  cout << " Look for rising and folling edges of logic time box: " << endl;
178  int start = -1;
179  int lenght = -1;
180  for(int j = 1; j < nBins;j++) {
181  int diff = (int)hLTB.GetBinContent(j+1)-(int)hLTB.GetBinContent(j);
182  if(diff == 1) { // This is a rising edge
183  start = j;
184  lenght = 1;
185  if(theVerbosityLevel >= 2) {
186  cout << " >>>----" << endl;
187  cout << " bin: " << j << " is a rising edge" << endl;
188  }
189  } else if(diff == -1) { // This is a falling edge
190  startAndLenght.push_back(make_pair(start, lenght));
191  if(theVerbosityLevel >= 2) {
192  cout << " bin: " << j << " is a falling edge, lenght is: " << lenght << endl;
193  cout << " <<<----" << endl;
194  }
195  start = -1;
196  lenght = -1;
197  } else if(diff == 0 && (int)hLTB.GetBinContent(j) == 1) { // This is a bin within the plateau
198  lenght ++;
199  }
200  }
201 
202  if(theVerbosityLevel >= 2)
203  cout << " Add macro intervals: " << endl;
204  // Look for consecutive plateau: 2 plateau are consecutive if the gap is < 5 ns
205  vector< pair<int, int> > superIntervals;
206  for(vector< pair<int, int> >::const_iterator interval = startAndLenght.begin();
207  interval != startAndLenght.end();
208  ++interval) {
209  pair<int, int> theInterval = (*interval);
210  vector< pair<int, int> >::const_iterator next = interval;
211  while(++next != startAndLenght.end()) {
212  int gap = (*next).first - (theInterval.first+theInterval.second);
213  double gabInNs = binValue*gap;
214  if(theVerbosityLevel >= 2)
215  cout << " gap: " << gabInNs << "(ns)" << endl;
216  if(gabInNs > 20) {
217  break;
218  } else {
219  theInterval = make_pair(theInterval.first, theInterval.second+gap+(*next).second);
220  superIntervals.push_back(theInterval);
221  if(theVerbosityLevel >= 2)
222  cout << " Add interval, start: " << theInterval.first
223  << " width: " << theInterval.second << endl;
224  }
225  }
226  }
227 
228  // merge the vectors of intervals
229  copy(superIntervals.begin(), superIntervals.end(), back_inserter(startAndLenght));
230 
231  // Look for the plateau of the right lenght
232  if(theVerbosityLevel >= 2)
233  cout << " Look for the best interval:" << endl;
234  int delta = 999999;
235  int beginning = -1;
236  int tbWidth = -1;
237  for(vector< pair<int, int> >::const_iterator stAndL = startAndLenght.begin();
238  stAndL != startAndLenght.end();
239  ++stAndL) {
240  if(abs((*stAndL).second - tBoxWidth) < delta) {
241  delta = abs((*stAndL).second - tBoxWidth);
242  beginning = (*stAndL).first;
243  tbWidth = (*stAndL).second;
244  if(theVerbosityLevel >= 2)
245  cout << " Candidate: First: " << beginning
246  << ", width: " << tbWidth
247  << ", delta: " << delta << endl;
248  }
249  }
250 
251  mean = xMin + beginning*binValue;
252  sigma = theSigma; //FIXME: estimate it!
253 
254  tBoxMax = hTBox->GetMaximum();
255 
256  // Define the fit range
257  xFitMin = mean-5.*sigma;
258  xFitMax = mean+5.*sigma;
259 
260  if(theVerbosityLevel >= 1) {
261  cout << " Time Box Rising edge: " << mean << endl;
262  cout << " Time Box Width: " << tbWidth*binValue << endl;
263  cout << " = Seeds and range for fit:" << endl;
264  cout << " Seed mean = " << mean << endl;
265  cout << " Seed sigma = " << sigma << endl;
266  cout << " Fitting from = " << xFitMin << " to " << xFitMax << endl << endl;
267  }
268 }
Definition: start.py:1
dbl * delta
Definition: mlp_gen.cc:36
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
unsigned int theVerbosityLevel
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 107 of file DTTimeBoxFitter.cc.

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

Referenced by fitTimeBox().

108  {
109  if(theVerbosityLevel >= 1)
110  cout << " === Insert seeds for the Time Box fit:" << endl;
111 
112  cout << "Inser the fit mean:" << endl;
113  cin >> mean;
114 
115  sigma = theSigma; //FIXME: estimate it!
116 
117  tBoxMax = hTBox->GetMaximum();
118 
119  // Define the fit range
120  xFitMin = mean-5.*sigma;
121  xFitMax = mean+5.*sigma;
122 
123  if(theVerbosityLevel >= 1) {
124  cout << " Time Box Rising edge: " << mean << endl;
125  cout << " = Seeds and range for fit:" << endl;
126  cout << " Seed mean = " << mean << endl;
127  cout << " Seed sigma = " << sigma << endl;
128  cout << " Fitting from = " << xFitMin << " to " << xFitMax << endl << endl;
129  }
130 }
unsigned int theVerbosityLevel
void DTTimeBoxFitter::setFitSigma ( double  sigma)
inline

Definition at line 56 of file DTTimeBoxFitter.h.

References theSigma.

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

Switch to interactive fit.

Definition at line 47 of file DTTimeBoxFitter.h.

References interactiveFit.

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

Set the rebin.

Definition at line 52 of file DTTimeBoxFitter.h.

References rebin.

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

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

Definition at line 42 of file DTTimeBoxFitter.h.

References theVerbosityLevel.

42  {
43  theVerbosityLevel = lvl;
44  }
unsigned int theVerbosityLevel

Member Data Documentation

TFile* DTTimeBoxFitter::hDebugFile
private

Definition at line 64 of file DTTimeBoxFitter.h.

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

bool DTTimeBoxFitter::interactiveFit
private

Definition at line 67 of file DTTimeBoxFitter.h.

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

int DTTimeBoxFitter::rebin
private

Definition at line 68 of file DTTimeBoxFitter.h.

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

double DTTimeBoxFitter::theSigma
private

Definition at line 69 of file DTTimeBoxFitter.h.

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

unsigned int DTTimeBoxFitter::theVerbosityLevel
private

Definition at line 66 of file DTTimeBoxFitter.h.

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