27 if(debugFileName !=
"")
hDebugFile =
new TFile(debugFileName.Data(),
"RECREATE");
44 hTimeBox->Rebin(
rebin);
47 if(hTimeBox->GetEntries() == 0) {
48 cout <<
"[DTTimeBoxFitter]***Error: the time box contains no entry!" << endl;
49 return make_pair(-1, -1);
52 if(hTimeBox->GetEntries() < 50000) {
66 TH1F *hTimeBoxForSeed = (TH1F*) hTimeBox->Clone();
68 getFitSeeds(hTimeBoxForSeed, xValue, xFitSigma, tBoxMax, xFitMin, xFitMax);
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);
89 hTimeBox->Fit(
"IntGauss", option.c_str(),
"",xFitMin, xFitMax);
92 double mean = fIntGaus->GetParameter(
"Mean");
93 double sigma = fIntGaus->GetParameter(
"Sigma");
95 double chiSquare = fIntGaus->GetChisquare()/fIntGaus->GetNDF();
98 cout <<
" === Fit Results: " << endl;
99 cout <<
" Fitted mean = " << mean << endl;
100 cout <<
" Fitted sigma = " << sigma << endl;
101 cout <<
" Reduced Chi Square = " << chiSquare << endl;
103 return make_pair(mean, sigma);
110 double& xFitMin,
double& xFitMax) {
112 cout <<
" === Insert seeds for the Time Box fit:" << endl;
114 cout <<
"Inser the fit mean:" << endl;
119 tBoxMax = hTBox->GetMaximum();
122 xFitMin = mean-5.*sigma;
123 xFitMax = mean+5.*sigma;
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;
137 double& xFitMin,
double& xFitMax) {
139 cout <<
" === Looking for fit seeds in Time Box:" << endl;
143 static const int tBoxWidth = 400;
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();
150 double binValue = (double)(xMax-xMin)/(double)nBins;
153 const double threshold = binValue*nEntries/(double)(tBoxWidth*3.);
155 cout <<
" Threshold for logic time box is (# entries): " << threshold << endl;
158 while(threshold > hTBox->GetMaximum()/2. && nRebins < 5) {
159 cout <<
" Rebinning!" << endl;
161 nBins = hTBox->GetNbinsX();
162 binValue = (double)(xMax-xMin)/(double)nBins;
167 TString hLName = TString(hTBox->GetName())+
"L";
168 TH1F hLTB(hLName.Data(),
"Logic Time Box", nBins, xMin, xMax);
170 for(
int i = 1;
i <= nBins;
i++) {
172 hLTB.SetBinContent(
i, 1);
177 vector< pair<int, int> > startAndLenght;
179 cout <<
" Look for rising and folling edges of logic time box: " << endl;
182 for(
int j = 1;
j < nBins;
j++) {
183 int diff = (int)hLTB.GetBinContent(
j+1)-(int)hLTB.GetBinContent(
j);
188 cout <<
" >>>----" << endl;
189 cout <<
" bin: " <<
j <<
" is a rising edge" << endl;
191 }
else if(diff == -1) {
192 startAndLenght.push_back(make_pair(start, lenght));
194 cout <<
" bin: " <<
j <<
" is a falling edge, lenght is: " << lenght << endl;
195 cout <<
" <<<----" << endl;
199 }
else if(diff == 0 && (
int)hLTB.GetBinContent(
j) == 1) {
205 cout <<
" Add macro intervals: " << endl;
207 vector< pair<int, int> > superIntervals;
208 for(vector< pair<int, int> >::const_iterator
interval = startAndLenght.begin();
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;
217 cout <<
" gap: " << gabInNs <<
"(ns)" << endl;
221 theInterval = make_pair(theInterval.first, theInterval.second+gap+(*next).second);
222 superIntervals.push_back(theInterval);
224 cout <<
" Add interval, start: " << theInterval.first
225 <<
" width: " << theInterval.second << endl;
231 copy(superIntervals.begin(), superIntervals.end(), back_inserter(startAndLenght));
235 cout <<
" Look for the best interval:" << endl;
239 for(vector< pair<int, int> >::const_iterator stAndL = startAndLenght.begin();
240 stAndL != startAndLenght.end();
242 if(
abs((*stAndL).second - tBoxWidth) <
delta) {
243 delta =
abs((*stAndL).second - tBoxWidth);
244 beginning = (*stAndL).first;
245 tbWidth = (*stAndL).second;
247 cout <<
" Candidate: First: " << beginning
248 <<
", width: " << tbWidth
249 <<
", delta: " << delta << endl;
253 mean = xMin + beginning*binValue;
256 tBoxMax = hTBox->GetMaximum();
259 xFitMin = mean-5.*sigma;
260 xFitMax = mean+5.*sigma;
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;
275 double media = par[1];
276 double sigma = par[2];
277 double var = (x[0]-media)/(sigma*
sqrt(2.));
279 return 0.5*par[0]*(1+TMath::Erf(var));
tuple start
Check for commandline option errors.
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.
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) ...