File indexing completed on 2024-04-06 11:58:30
0001
0002
0003
0004
0005
0006
0007
0008 #include "CalibMuon/DTCalibration/interface/DTTimeBoxFitter.h"
0009
0010 #include <iostream>
0011 #include <vector>
0012
0013 #include "TFile.h"
0014 #include "TH1F.h"
0015 #include "TMath.h"
0016 #include "TF1.h"
0017 #include "TString.h"
0018
0019 using namespace std;
0020
0021 DTTimeBoxFitter::DTTimeBoxFitter(const TString& debugFileName)
0022 : hDebugFile(nullptr), theVerbosityLevel(0), theSigma(10.) {
0023
0024 if (debugFileName != "")
0025 hDebugFile = new TFile(debugFileName.Data(), "RECREATE");
0026 interactiveFit = false;
0027 rebin = 1;
0028 }
0029
0030 DTTimeBoxFitter::~DTTimeBoxFitter() {
0031 if (hDebugFile != nullptr)
0032 hDebugFile->Close();
0033 }
0034
0035
0036 pair<double, double> DTTimeBoxFitter::fitTimeBox(TH1F* hTimeBox) {
0037 hTimeBox->Rebin(rebin);
0038
0039
0040 if (hTimeBox->GetEntries() == 0) {
0041 cout << "[DTTimeBoxFitter]***Error: the time box contains no entry!" << endl;
0042 return make_pair(-1, -1);
0043 }
0044
0045 if (hTimeBox->GetEntries() < 50000) {
0046 hTimeBox->Rebin(2);
0047 }
0048
0049
0050
0051 double xFitMin = 0;
0052 double xFitMax = 0;
0053 double xValue = 0;
0054 double xFitSigma = 0;
0055 double tBoxMax = 0;
0056
0057
0058
0059 TH1F* hTimeBoxForSeed = (TH1F*)hTimeBox->Clone();
0060 if (!interactiveFit) {
0061 getFitSeeds(hTimeBoxForSeed, xValue, xFitSigma, tBoxMax, xFitMin, xFitMax);
0062 } else {
0063 getInteractiveFitSeeds(hTimeBoxForSeed, xValue, xFitSigma, tBoxMax, xFitMin, xFitMax);
0064 }
0065
0066
0067 TF1* fIntGaus = new TF1("IntGauss", intGauss, xFitMin, xFitMax, 3);
0068 fIntGaus->SetParName(0, "Constant");
0069 fIntGaus->SetParameter(0, tBoxMax);
0070 fIntGaus->SetParName(1, "Mean");
0071 fIntGaus->SetParameter(1, xValue);
0072 fIntGaus->SetParName(2, "Sigma");
0073 fIntGaus->SetParameter(2, xFitSigma);
0074 fIntGaus->SetLineColor(kRed);
0075
0076
0077 string option = "Q";
0078 if (theVerbosityLevel >= 2)
0079 option = "";
0080
0081 hTimeBox->Fit(fIntGaus, option.c_str(), "", xFitMin, xFitMax);
0082
0083
0084 double mean = fIntGaus->GetParameter("Mean");
0085 double sigma = fIntGaus->GetParameter("Sigma");
0086
0087 double chiSquare = fIntGaus->GetChisquare() / fIntGaus->GetNDF();
0088
0089 if (theVerbosityLevel >= 1) {
0090 cout << " === Fit Results: " << endl;
0091 cout << " Fitted mean = " << mean << endl;
0092 cout << " Fitted sigma = " << sigma << endl;
0093 cout << " Reduced Chi Square = " << chiSquare << endl;
0094 }
0095 return make_pair(mean, sigma);
0096 }
0097
0098
0099 void DTTimeBoxFitter::getInteractiveFitSeeds(
0100 TH1F* hTBox, double& mean, double& sigma, double& tBoxMax, double& xFitMin, double& xFitMax) {
0101 if (theVerbosityLevel >= 1)
0102 cout << " === Insert seeds for the Time Box fit:" << endl;
0103
0104 cout << "Inser the fit mean:" << endl;
0105 cin >> mean;
0106
0107 sigma = theSigma;
0108
0109 tBoxMax = hTBox->GetMaximum();
0110
0111
0112 xFitMin = mean - 5. * sigma;
0113 xFitMax = mean + 5. * sigma;
0114
0115 if (theVerbosityLevel >= 1) {
0116 cout << " Time Box Rising edge: " << mean << endl;
0117 cout << " = Seeds and range for fit:" << endl;
0118 cout << " Seed mean = " << mean << endl;
0119 cout << " Seed sigma = " << sigma << endl;
0120 cout << " Fitting from = " << xFitMin << " to " << xFitMax << endl << endl;
0121 }
0122 }
0123
0124
0125 void DTTimeBoxFitter::getFitSeeds(
0126 TH1F* hTBox, double& mean, double& sigma, double& tBoxMax, double& xFitMin, double& xFitMax) {
0127 if (theVerbosityLevel >= 1)
0128 cout << " === Looking for fit seeds in Time Box:" << endl;
0129
0130
0131 static const int tBoxWidth = 400;
0132
0133 int nBins = hTBox->GetNbinsX();
0134 const int xMin = (int)hTBox->GetXaxis()->GetXmin();
0135 const int xMax = (int)hTBox->GetXaxis()->GetXmax();
0136 const int nEntries = (int)hTBox->GetEntries();
0137
0138 double binValue = (double)(xMax - xMin) / (double)nBins;
0139
0140
0141 const double threshold = binValue * nEntries / (double)(tBoxWidth * 3.);
0142 if (theVerbosityLevel >= 2)
0143 cout << " Threshold for logic time box is (# entries): " << threshold << endl;
0144
0145 int nRebins = 0;
0146 while (threshold > hTBox->GetMaximum() / 2. && nRebins < 5) {
0147 cout << " Rebinning!" << endl;
0148 hTBox->Rebin(2);
0149 nBins = hTBox->GetNbinsX();
0150 binValue = (double)(xMax - xMin) / (double)nBins;
0151 nRebins++;
0152 }
0153
0154 if (hDebugFile != nullptr)
0155 hDebugFile->cd();
0156 TString hLName = TString(hTBox->GetName()) + "L";
0157 TH1F hLTB(hLName.Data(), "Logic Time Box", nBins, xMin, xMax);
0158
0159 for (int i = 1; i <= nBins; i++) {
0160 if (hTBox->GetBinContent(i) > threshold)
0161 hLTB.SetBinContent(i, 1);
0162 }
0163 if (hDebugFile != nullptr)
0164 hLTB.Write();
0165
0166
0167 vector<pair<int, int> > startAndLenght;
0168 if (theVerbosityLevel >= 2)
0169 cout << " Look for rising and folling edges of logic time box: " << endl;
0170 int start = -1;
0171 int lenght = -1;
0172 for (int j = 1; j < nBins; j++) {
0173 int diff = (int)hLTB.GetBinContent(j + 1) - (int)hLTB.GetBinContent(j);
0174 if (diff == 1) {
0175 start = j;
0176 lenght = 1;
0177 if (theVerbosityLevel >= 2) {
0178 cout << " >>>----" << endl;
0179 cout << " bin: " << j << " is a rising edge" << endl;
0180 }
0181 } else if (diff == -1) {
0182 startAndLenght.push_back(make_pair(start, lenght));
0183 if (theVerbosityLevel >= 2) {
0184 cout << " bin: " << j << " is a falling edge, lenght is: " << lenght << endl;
0185 cout << " <<<----" << endl;
0186 }
0187 start = -1;
0188 lenght = -1;
0189 } else if (diff == 0 && (int)hLTB.GetBinContent(j) == 1) {
0190 lenght++;
0191 }
0192 }
0193
0194 if (theVerbosityLevel >= 2)
0195 cout << " Add macro intervals: " << endl;
0196
0197 vector<pair<int, int> > superIntervals;
0198 for (vector<pair<int, int> >::const_iterator interval = startAndLenght.begin(); interval != startAndLenght.end();
0199 ++interval) {
0200 pair<int, int> theInterval = (*interval);
0201 vector<pair<int, int> >::const_iterator next = interval;
0202 while (++next != startAndLenght.end()) {
0203 int gap = (*next).first - (theInterval.first + theInterval.second);
0204 double gabInNs = binValue * gap;
0205 if (theVerbosityLevel >= 2)
0206 cout << " gap: " << gabInNs << "(ns)" << endl;
0207 if (gabInNs > 20) {
0208 break;
0209 } else {
0210 theInterval = make_pair(theInterval.first, theInterval.second + gap + (*next).second);
0211 superIntervals.push_back(theInterval);
0212 if (theVerbosityLevel >= 2)
0213 cout << " Add interval, start: " << theInterval.first << " width: " << theInterval.second << endl;
0214 }
0215 }
0216 }
0217
0218
0219 copy(superIntervals.begin(), superIntervals.end(), back_inserter(startAndLenght));
0220
0221
0222 if (theVerbosityLevel >= 2)
0223 cout << " Look for the best interval:" << endl;
0224 int delta = 999999;
0225 int beginning = -1;
0226 int tbWidth = -1;
0227 for (vector<pair<int, int> >::const_iterator stAndL = startAndLenght.begin(); stAndL != startAndLenght.end();
0228 ++stAndL) {
0229 if (abs((*stAndL).second - tBoxWidth) < delta) {
0230 delta = abs((*stAndL).second - tBoxWidth);
0231 beginning = (*stAndL).first;
0232 tbWidth = (*stAndL).second;
0233 if (theVerbosityLevel >= 2)
0234 cout << " Candidate: First: " << beginning << ", width: " << tbWidth << ", delta: " << delta << endl;
0235 }
0236 }
0237
0238 mean = xMin + beginning * binValue;
0239 sigma = theSigma;
0240
0241 tBoxMax = hTBox->GetMaximum();
0242
0243
0244 xFitMin = mean - 5. * sigma;
0245 xFitMax = mean + 5. * sigma;
0246
0247 if (theVerbosityLevel >= 1) {
0248 cout << " Time Box Rising edge: " << mean << endl;
0249 cout << " Time Box Width: " << tbWidth * binValue << endl;
0250 cout << " = Seeds and range for fit:" << endl;
0251 cout << " Seed mean = " << mean << endl;
0252 cout << " Seed sigma = " << sigma << endl;
0253 cout << " Fitting from = " << xFitMin << " to " << xFitMax << endl << endl;
0254 }
0255 }
0256
0257 double intGauss(double* x, double* par) {
0258 double media = par[1];
0259 double sigma = par[2];
0260 double var = (x[0] - media) / (sigma * sqrt(2.));
0261
0262 return 0.5 * par[0] * (1 + TMath::Erf(var));
0263 }