#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 17 of file DTTimeBoxFitter.h.

Constructor & Destructor Documentation

◆ DTTimeBoxFitter()

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

Constructor.

Definition at line 21 of file DTTimeBoxFitter.cc.

References hDebugFile, interactiveFit, and rebin.

     : hDebugFile(nullptr), theVerbosityLevel(0), theSigma(10.) {
   // Create a root file for debug output only if needed
   if (debugFileName != "")
     hDebugFile = new TFile(debugFileName.Data(), "RECREATE");
   interactiveFit = false;
   rebin = 1;
 }

◆ ~DTTimeBoxFitter()

DTTimeBoxFitter::~DTTimeBoxFitter ( )

virtual

Destructor.

Definition at line 30 of file DTTimeBoxFitter.cc.

References hDebugFile.

                                   {
   if (hDebugFile != nullptr)
     hDebugFile->Close();
 }

Member Function Documentation

◆ fitTimeBox()

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

References simKBmtfDigis_cfi::chiSquare, gather_cfg::cout, getFitSeeds(), getInteractiveFitSeeds(), interactiveFit, intGauss(), SiStripPI::mean, fileinputsource_cfi::option, rebin, and theVerbosityLevel.

Referenced by DTTTrigWriter::analyze().

                                                                {
   hTimeBox->Rebin(rebin);
 
   // Check if the histo contains any entry
   if (hTimeBox->GetEntries() == 0) {
     cout << "[DTTimeBoxFitter]***Error: the time box contains no entry!" << endl;
     return make_pair(-1, -1);
   }
 
   if (hTimeBox->GetEntries() < 50000) {
     hTimeBox->Rebin(2);
   }
 
   // Get seeds for the fit
   // The TimeBox range to be fitted (the rising edge)
   double xFitMin = 0;    // Min value for the fit
   double xFitMax = 0;    // Max value for the fit
   double xValue = 0;     // The mean value of the gaussian
   double xFitSigma = 0;  // The sigma of the gaussian
   double tBoxMax = 0;    // The max of the time box, it is used as seed for gaussian integral
 
   //hTimeBox->Rebin(2); //FIXME: Temporary for low statistics
 
   TH1F* hTimeBoxForSeed = (TH1F*)hTimeBox->Clone();  //FIXME: test
   if (!interactiveFit) {
     getFitSeeds(hTimeBoxForSeed, xValue, xFitSigma, tBoxMax, xFitMin, xFitMax);
   } else {
     getInteractiveFitSeeds(hTimeBoxForSeed, xValue, xFitSigma, tBoxMax, xFitMin, xFitMax);
   }
 
   // Define the fitting function and use fit seeds
   TF1* fIntGaus = new TF1("IntGauss", intGauss, xFitMin, xFitMax, 3);
   fIntGaus->SetParName(0, "Constant");
   fIntGaus->SetParameter(0, tBoxMax);
   fIntGaus->SetParName(1, "Mean");
   fIntGaus->SetParameter(1, xValue);
   fIntGaus->SetParName(2, "Sigma");
   fIntGaus->SetParameter(2, xFitSigma);
   fIntGaus->SetLineColor(kRed);
 
   // Fit the histo
   string option = "Q";
   if (theVerbosityLevel >= 2)
     option = "";
 
   hTimeBox->Fit(fIntGaus, option.c_str(), "", xFitMin, xFitMax);
 
   // Get fitted parameters
   double mean = fIntGaus->GetParameter("Mean");
   double sigma = fIntGaus->GetParameter("Sigma");
   //   double constant = fIntGaus->GetParameter("Constant");
   double chiSquare = fIntGaus->GetChisquare() / fIntGaus->GetNDF();
 
   if (theVerbosityLevel >= 1) {
     cout << " === Fit Results: " << endl;
     cout << "     Fitted mean = " << mean << endl;
     cout << "     Fitted sigma = " << sigma << endl;
     cout << "     Reduced Chi Square = " << chiSquare << endl;
   }
   return make_pair(mean, sigma);
 }

◆ getFitSeeds()

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

References funct::abs(), filterCSVwithJSON::copy, gather_cfg::cout, dumpMFGeometry_cfg::delta, change_name::diff, HLT_2022v15_cff::gap, hDebugFile, mps_fire::i, createfilelist::int, writeEcalDQMStatus::interval, dqmiolumiharvest::j, SiStripPI::mean, seedmultiplicitymonitor_newtracking_cfi::nBins, GetRecoTauVFromDQM_MC_cff::next, theSigma, theVerbosityLevel, remoteMonitoring_LASER_era2018_cfg::threshold, trackerHitRTTI::vector, multiplicitycorr_cfi::xMax, and photonAnalyzer_cfi::xMin.

Referenced by fitTimeBox().

                                                                                                  {
   if (theVerbosityLevel >= 1)
     cout << " === Looking for fit seeds in Time Box:" << endl;
 
   // The approximate width of the time box
   static const int tBoxWidth = 400;  //FIXE: tune it
 
   int nBins = hTBox->GetNbinsX();
   const int xMin = (int)hTBox->GetXaxis()->GetXmin();
   const int xMax = (int)hTBox->GetXaxis()->GetXmax();
   const int nEntries = (int)hTBox->GetEntries();
 
   double binValue = (double)(xMax - xMin) / (double)nBins;
 
   // Compute a threshold for TimeBox discrimination
   const double threshold = binValue * nEntries / (double)(tBoxWidth * 3.);
   if (theVerbosityLevel >= 2)
     cout << "   Threshold for logic time box is (# entries): " << threshold << endl;
 
   int nRebins = 0;  // protection for infinite loop
   while (threshold > hTBox->GetMaximum() / 2. && nRebins < 5) {
     cout << " Rebinning!" << endl;
     hTBox->Rebin(2);
     nBins = hTBox->GetNbinsX();
     binValue = (double)(xMax - xMin) / (double)nBins;
     nRebins++;
   }
 
   if (hDebugFile != nullptr)
     hDebugFile->cd();
   TString hLName = TString(hTBox->GetName()) + "L";
   TH1F hLTB(hLName.Data(), "Logic Time Box", nBins, xMin, xMax);
   // Loop over all time box bins and discriminate them accordigly to the threshold
   for (int i = 1; i <= nBins; i++) {
     if (hTBox->GetBinContent(i) > threshold)
       hLTB.SetBinContent(i, 1);
   }
   if (hDebugFile != nullptr)
     hLTB.Write();
 
   // Look for the time box in the "logic histo" and save beginning and lenght of each plateau
   vector<pair<int, int> > startAndLenght;
   if (theVerbosityLevel >= 2)
     cout << "   Look for rising and folling edges of logic time box: " << endl;
   int start = -1;
   int lenght = -1;
   for (int j = 1; j < nBins; j++) {
     int diff = (int)hLTB.GetBinContent(j + 1) - (int)hLTB.GetBinContent(j);
     if (diff == 1) {  // This is a rising edge
       start = j;
       lenght = 1;
       if (theVerbosityLevel >= 2) {
         cout << "     >>>----" << endl;
         cout << "      bin: " << j << " is a rising edge" << endl;
       }
     } else if (diff == -1) {  // This is a falling edge
       startAndLenght.push_back(make_pair(start, lenght));
       if (theVerbosityLevel >= 2) {
         cout << "      bin: " << j << " is a falling edge, lenght is: " << lenght << endl;
         cout << "     <<<----" << endl;
       }
       start = -1;
       lenght = -1;
     } else if (diff == 0 && (int)hLTB.GetBinContent(j) == 1) {  // This is a bin within the plateau
       lenght++;
     }
   }
 
   if (theVerbosityLevel >= 2)
     cout << "   Add macro intervals: " << endl;
   // Look for consecutive plateau: 2 plateau are consecutive if the gap is < 5 ns
   vector<pair<int, int> > superIntervals;
   for (vector<pair<int, int> >::const_iterator interval = startAndLenght.begin(); interval != startAndLenght.end();
        ++interval) {
     pair<int, int> theInterval = (*interval);
     vector<pair<int, int> >::const_iterator next = interval;
     while (++next != startAndLenght.end()) {
       int gap = (*next).first - (theInterval.first + theInterval.second);
       double gabInNs = binValue * gap;
       if (theVerbosityLevel >= 2)
         cout << "      gap: " << gabInNs << "(ns)" << endl;
       if (gabInNs > 20) {
         break;
       } else {
         theInterval = make_pair(theInterval.first, theInterval.second + gap + (*next).second);
         superIntervals.push_back(theInterval);
         if (theVerbosityLevel >= 2)
           cout << "          Add interval, start: " << theInterval.first << " width: " << theInterval.second << endl;
       }
     }
   }
 
   // merge the vectors of intervals
   copy(superIntervals.begin(), superIntervals.end(), back_inserter(startAndLenght));
 
   // Look for the plateau of the right lenght
   if (theVerbosityLevel >= 2)
     cout << "    Look for the best interval:" << endl;
   int delta = 999999;
   int beginning = -1;
   int tbWidth = -1;
   for (vector<pair<int, int> >::const_iterator stAndL = startAndLenght.begin(); stAndL != startAndLenght.end();
        ++stAndL) {
     if (abs((*stAndL).second - tBoxWidth) < delta) {
       delta = abs((*stAndL).second - tBoxWidth);
       beginning = (*stAndL).first;
       tbWidth = (*stAndL).second;
       if (theVerbosityLevel >= 2)
         cout << "   Candidate: First: " << beginning << ", width: " << tbWidth << ", delta: " << delta << endl;
     }
   }
 
   mean = xMin + beginning * binValue;
   sigma = theSigma;  //FIXME: estimate it!
 
   tBoxMax = hTBox->GetMaximum();
 
   // Define the fit range
   xFitMin = mean - 5. * sigma;
   xFitMax = mean + 5. * sigma;
 
   if (theVerbosityLevel >= 1) {
     cout << "      Time Box Rising edge: " << mean << endl;
     cout << "      Time Box Width: " << tbWidth * binValue << endl;
     cout << "    = Seeds and range for fit:" << endl;
     cout << "       Seed mean = " << mean << endl;
     cout << "       Seed sigma = " << sigma << endl;
     cout << "       Fitting from = " << xFitMin << " to " << xFitMax << endl << endl;
   }
 }

◆ getInteractiveFitSeeds()

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

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

Referenced by fitTimeBox().

                                                                                                  {
   if (theVerbosityLevel >= 1)
     cout << " === Insert seeds for the Time Box fit:" << endl;
 
   cout << "Inser the fit mean:" << endl;
   cin >> mean;
 
   sigma = theSigma;  //FIXME: estimate it!
 
   tBoxMax = hTBox->GetMaximum();
 
   // Define the fit range
   xFitMin = mean - 5. * sigma;
   xFitMax = mean + 5. * sigma;
 
   if (theVerbosityLevel >= 1) {
     cout << "      Time Box Rising edge: " << mean << endl;
     cout << "    = Seeds and range for fit:" << endl;
     cout << "       Seed mean = " << mean << endl;
     cout << "       Seed sigma = " << sigma << endl;
     cout << "       Fitting from = " << xFitMin << " to " << xFitMax << endl << endl;
   }
 }