DTCalibration/src/DTTimeBoxFitter.cc

0001
0002 /*
0003  *  See header file for a description of this class.
0004  *
0005  *  \author G. Cerminara - INFN Torino
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   // Create a root file for debug output only if needed
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 /// Compute the ttrig (in ns) from the Time Box
0036 pair<double, double> DTTimeBoxFitter::fitTimeBox(TH1F* hTimeBox) {
0037   hTimeBox->Rebin(rebin);
0038
0039   // Check if the histo contains any entry
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   // Get seeds for the fit
0050   // The TimeBox range to be fitted (the rising edge)
0051   double xFitMin = 0;    // Min value for the fit
0052   double xFitMax = 0;    // Max value for the fit
0053   double xValue = 0;     // The mean value of the gaussian
0054   double xFitSigma = 0;  // The sigma of the gaussian
0055   double tBoxMax = 0;    // The max of the time box, it is used as seed for gaussian integral
0056
0057   //hTimeBox->Rebin(2); //FIXME: Temporary for low statistics
0058
0059   TH1F* hTimeBoxForSeed = (TH1F*)hTimeBox->Clone();  //FIXME: test
0060   if (!interactiveFit) {
0061     getFitSeeds(hTimeBoxForSeed, xValue, xFitSigma, tBoxMax, xFitMin, xFitMax);
0062   } else {
0063     getInteractiveFitSeeds(hTimeBoxForSeed, xValue, xFitSigma, tBoxMax, xFitMin, xFitMax);
0064   }
0065
0066   // Define the fitting function and use fit seeds
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   // Fit the histo
0077   string option = "Q";
0078   if (theVerbosityLevel >= 2)
0079     option = "";
0080
0081   hTimeBox->Fit(fIntGaus, option.c_str(), "", xFitMin, xFitMax);
0082
0083   // Get fitted parameters
0084   double mean = fIntGaus->GetParameter("Mean");
0085   double sigma = fIntGaus->GetParameter("Sigma");
0086   //   double constant = fIntGaus->GetParameter("Constant");
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 // Get the seeds for the fit as input from user!It is used if interactiveFit == true
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;  //FIXME: estimate it!
0108
0109   tBoxMax = hTBox->GetMaximum();
0110
0111   // Define the fit range
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 // Automatically compute the seeds the range to be used for time box fit
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   // The approximate width of the time box
0131   static const int tBoxWidth = 400;  //FIXE: tune it
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   // Compute a threshold for TimeBox discrimination
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;  // protection for infinite loop
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   // Loop over all time box bins and discriminate them accordigly to the threshold
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   // Look for the time box in the "logic histo" and save beginning and lenght of each plateau
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) {  // This is a rising edge
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) {  // This is a falling edge
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) {  // This is a bin within the plateau
0190       lenght++;
0191     }
0192   }
0193
0194   if (theVerbosityLevel >= 2)
0195     cout << "   Add macro intervals: " << endl;
0196   // Look for consecutive plateau: 2 plateau are consecutive if the gap is < 5 ns
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   // merge the vectors of intervals
0219   copy(superIntervals.begin(), superIntervals.end(), back_inserter(startAndLenght));
0220
0221   // Look for the plateau of the right lenght
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;  //FIXME: estimate it!
0240
0241   tBoxMax = hTBox->GetMaximum();
0242
0243   // Define the fit range
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 }