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.

                                                             : hDebugFile(0),
                                 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 ( )

virtual

Destructor.

Definition at line 32 of file DTTimeBoxFitter.cc.

References hDebugFile.

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

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, ctppsDiamondLocalTracks_cfi::sigma, and theVerbosityLevel.

                                                               {

    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);
}

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, ctppsDiamondLocalTracks_cfi::sigma, theSigma, theVerbosityLevel, ctppsDiamondLocalTracks_cfi::threshold, anotherprimaryvertexanalyzer_cfi::xMax, and anotherprimaryvertexanalyzer_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 != 0) 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 != 0) 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;
  }
}

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, ctppsDiamondLocalTracks_cfi::sigma, 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;
  }
}

void DTTimeBoxFitter::setFitSigma ( double  sigma )

inline

Definition at line 56 of file DTTimeBoxFitter.h.

References ctppsDiamondLocalTracks_cfi::sigma, and theSigma.

                                 {
    theSigma = sigma;
  }

void DTTimeBoxFitter::setInteractiveFit ( bool  isInteractive )

inline

Switch to interactive fit.

Definition at line 47 of file DTTimeBoxFitter.h.

References interactiveFit.

                                             {
    interactiveFit = isInteractive;
  }

void DTTimeBoxFitter::setRebinning ( int  reb )

inline

Set the rebin.

Definition at line 52 of file DTTimeBoxFitter.h.

References rebin.

                             {
    rebin = reb;
  }

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.

                                      {
    theVerbosityLevel = lvl;
  }

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().