DTMeanTimerFitter Class Reference

#include <DTMeanTimerFitter.h>

Public Member Functions
	DTMeanTimerFitter (TFile *file)
	Constructor. More...
std::vector< float >	evaluateVDriftAndReso (const TString &N)
	Fit the TMax histos and evaluate VDrift and resolution. More...
TF1 *	fitTMax (TH1F *histo)
	Really do the fit. More...
void	setVerbosity (unsigned int lvl)
	Set the verbosity of the output: 0 = silent, 1 = info, 2 = debug. More...
virtual	~DTMeanTimerFitter ()
	Destructor. More...

Private Attributes
TFile *	hDebugFile
TFile *	hInputFile
unsigned int	theVerbosityLevel

Detailed Description

Fit the Tmax histograms with a gaussian returning the mean values and the sigmas.

Author

S. Bolognesi - INFN Torino

Definition at line 17 of file DTMeanTimerFitter.h.

Constructor & Destructor Documentation

DTMeanTimerFitter()

DTMeanTimerFitter::DTMeanTimerFitter

(

TFile *

file

)

Constructor.

Definition at line 21 of file DTMeanTimerFitter.cc.

References hDebugFile.

                                                : hInputFile(file), theVerbosityLevel(0) {
  //hInputFile = new TFile("DTTMaxHistos.root", "READ");
  hDebugFile = new TFile("DTMeanTimerFitter.root", "RECREATE");
}

~DTMeanTimerFitter()

DTMeanTimerFitter::~DTMeanTimerFitter ( )

virtual

Destructor.

Definition at line 26 of file DTMeanTimerFitter.cc.

References hDebugFile.

{ hDebugFile->Close(); }

Member Function Documentation

evaluateVDriftAndReso()

vector< float > DTMeanTimerFitter::evaluateVDriftAndReso

(

const TString &

N

)

Fit the TMax histos and evaluate VDrift and resolution.

Definition at line 28 of file DTMeanTimerFitter.cc.

References submitPVResolutionJobs::count, cppFunctionSkipper::exception, validate-o2o-wbm::f1, fitTMax(), hDebugFile, hInputFile, mergeVDriftHistosByStation::histos, mps_fire::i, SiStripPI::mean, N, and mathSSE::sqrt().

                                                                       {
  // Retrieve histogram sets
  hTMaxCell *histos = new hTMaxCell(N, hInputFile);
  vector<float> vDriftAndReso;

  // Check that the histo for this cell exists
  if (histos->hTmax123 != nullptr) {
    vector<TH1F *> hTMax;  // histograms for <T_max> calculation
    vector<TH1F *> hT0;    // histograms for T0 evaluation
    hTMax.push_back(histos->hTmax123);
    hTMax.push_back(histos->hTmax124s72);
    hTMax.push_back(histos->hTmax124s78);
    hTMax.push_back(histos->hTmax134s72);
    hTMax.push_back(histos->hTmax134s78);
    hTMax.push_back(histos->hTmax234);

    hT0.push_back(histos->hTmax_3t0);
    hT0.push_back(histos->hTmax_2t0);
    hT0.push_back(histos->hTmax_t0);
    hT0.push_back(histos->hTmax_0);

    vector<Double_t> factor;          // factor relating the width of the Tmax distribution
                                      // and the cell resolution
    factor.push_back(sqrt(2. / 3.));  // hTmax123
    factor.push_back(sqrt(2. / 7.));  // hTmax124s72
    factor.push_back(sqrt(8. / 7.));  // hTmax124s78
    factor.push_back(sqrt(2. / 7.));  // hTmax134s72
    factor.push_back(sqrt(8. / 7.));  // hTmax134s78
    factor.push_back(sqrt(2. / 3.));  // hTmax234

    // Retrieve the gaussian mean and sigma for each TMax histogram
    vector<Double_t> mean;
    vector<Double_t> sigma;
    vector<Double_t> count;  //number of entries

    for (vector<TH1F *>::const_iterator ith = hTMax.begin(); ith != hTMax.end(); ++ith) {
      TF1 *funct = fitTMax(*ith);
      if (!funct) {
        edm::LogError("DTMeanTimerFitter") << "Error when fitting TMax..histogram name" << (*ith)->GetName();
        // return empty or -1 filled vector?
        vector<float> defvec(6, -1);
        return defvec;
      }
      hDebugFile->cd();
      (*ith)->Write();

      // Get mean, sigma and number of entries of each histogram
      mean.push_back(funct->GetParameter(1));
      sigma.push_back(funct->GetParameter(2));
      count.push_back((*ith)->GetEntries());
    }

    Double_t tMaxMean = 0.;
    Double_t wTMaxSum = 0.;
    Double_t sigmaT = 0.;
    Double_t wSigmaSum = 0.;

    //calculate total mean and sigma
    for (int i = 0; i <= 5; i++) {
      if (count[i] < 200)
        continue;
      tMaxMean += mean[i] * (count[i] / (sigma[i] * sigma[i]));
      wTMaxSum += count[i] / (sigma[i] * sigma[i]);
      sigmaT += count[i] * factor[i] * sigma[i];
      wSigmaSum += count[i];
      // cout << "TMaxMean "<<i<<": "<< mean[i] << " entries: " << count[i]
      // << " sigma: " << sigma[i]
      // << " weight: " << (count[i]/(sigma[i]*sigma[i])) << endl;
    }
    if ((!wTMaxSum) || (!wSigmaSum)) {
      edm::LogError("DTMeanTimerFitter") << "Error zero sum of weights..returning default";
      vector<float> defvec(6, -1);
      return defvec;
    }

    tMaxMean /= wTMaxSum;
    sigmaT /= wSigmaSum;

    //calculate v_drift and resolution
    Double_t vDrift = 2.1 / tMaxMean;  //2.1 is the half cell length in cm
    Double_t reso = vDrift * sigmaT;
    vDriftAndReso.push_back(vDrift);
    vDriftAndReso.push_back(reso);
    //if(theVerbosityLevel >= 1)
    edm::LogVerbatim("DTMeanTimerFitter")
        << " final TMaxMean=" << tMaxMean << " sigma= " << sigmaT << " v_d and reso: " << vDrift << " " << reso << endl;

    // Order t0 histogram by number of entries (choose histograms with higher nr. of entries)
    map<Double_t, TH1F *> hEntries;
    for (vector<TH1F *>::const_iterator ith = hT0.begin(); ith != hT0.end(); ++ith) {
      hEntries[(*ith)->GetEntries()] = (*ith);
    }

    // add at the end of hT0 the two hists with the higher number of entries
    int counter = 0;
    for (map<Double_t, TH1F *>::reverse_iterator iter = hEntries.rbegin(); iter != hEntries.rend(); ++iter) {
      counter++;
      if (counter == 1)
        hT0.push_back(iter->second);
      else if (counter == 2) {
        hT0.push_back(iter->second);
        break;
      }
    }

    // Retrieve the gaussian mean and sigma of histograms for Delta(t0) evaluation
    vector<Double_t> meanT0;
    vector<Double_t> sigmaT0;
    vector<Double_t> countT0;

    for (vector<TH1F *>::const_iterator ith = hT0.begin(); ith != hT0.end(); ++ith) {
      try {
        (*ith)->Fit("gaus");
      } catch (std::exception const &) {
        edm::LogError("DTMeanTimerFitter") << "Exception when fitting T0..histogram " << (*ith)->GetName();
        // return empty or -1 filled vector?
        vector<float> defvec(6, -1);
        return defvec;
      }
      TF1 *f1 = (*ith)->GetFunction("gaus");
      // Get mean, sigma and number of entries of the  histograms
      meanT0.push_back(f1->GetParameter(1));
      sigmaT0.push_back(f1->GetParameter(2));
      countT0.push_back((*ith)->GetEntries());
    }
    //calculate Delta(t0)
    if (hT0.size() != 6) {  // check if you have all the t0 hists
      edm::LogVerbatim("DTMeanTimerFitter") << "t0 histograms = " << hT0.size();
      for (int i = 1; i <= 4; i++) {
        vDriftAndReso.push_back(-1);
      }
      return vDriftAndReso;
    }

    for (int it0 = 0; it0 <= 2; it0++) {
      if ((countT0[it0] > 200) && (countT0[it0 + 1] > 200)) {
        Double_t deltaT0 = meanT0[it0] - meanT0[it0 + 1];
        vDriftAndReso.push_back(deltaT0);
      } else
        vDriftAndReso.push_back(999.);
    }
    //deltat0 using hists with max nr. of entries
    if ((countT0[4] > 200) && (countT0[5] > 200)) {
      Double_t t0Diff = histos->GetT0Factor(hT0[4]) - histos->GetT0Factor(hT0[5]);
      Double_t deltaT0MaxEntries = (meanT0[4] - meanT0[5]) / t0Diff;
      vDriftAndReso.push_back(deltaT0MaxEntries);
    } else
      vDriftAndReso.push_back(999.);
  } else {
    for (int i = 1; i <= 6; i++) {
      // 0=vdrift, 1=reso,  2=(3deltat0-2deltat0), 3=(2deltat0-1deltat0),
      // 4=(1deltat0-0deltat0), 5=deltat0 from hists with max entries,
      vDriftAndReso.push_back(-1);
    }
  }
  return vDriftAndReso;
}

fitTMax()

TF1 * DTMeanTimerFitter::fitTMax

(

TH1F *

histo

)

Really do the fit.

Definition at line 186 of file DTMeanTimerFitter.cc.

References cppFunctionSkipper::exception, validate-o2o-wbm::f1, timingPdfMaker::histo, LogDebug, and isotrackApplyRegressor::range.

Referenced by evaluateVDriftAndReso().

                                           {
  // Find distribution peak and fit range
  Double_t peak = (((((histo->GetXaxis())->GetXmax()) - ((histo->GetXaxis())->GetXmin())) / histo->GetNbinsX()) *
                   (histo->GetMaximumBin())) +
                  ((histo->GetXaxis())->GetXmin());
  //if(theVerbosityLevel >= 1)
  LogDebug("DTMeanTimerFitter") << "Peak " << peak << " : "
                                << "xmax " << ((histo->GetXaxis())->GetXmax()) << "            xmin "
                                << ((histo->GetXaxis())->GetXmin()) << "            nbin " << histo->GetNbinsX()
                                << "            bin with max " << (histo->GetMaximumBin());
  Double_t range = 2. * histo->GetRMS();

  // Fit each Tmax histogram with a Gaussian in a restricted interval
  TF1 *rGaus = new TF1("rGaus", "gaus", peak - range, peak + range);
  rGaus->SetMarkerSize();  // just silence gcc complainining about unused var
  try {
    histo->Fit("rGaus", "R");
  } catch (std::exception const &) {
    edm::LogError("DTMeanTimerFitter") << "Exception when fitting TMax..histogram " << histo->GetName()
                                       << "   setting return function pointer to zero";
    return nullptr;
  }
  TF1 *f1 = histo->GetFunction("rGaus");
  return f1;
}

setVerbosity()

void DTMeanTimerFitter::setVerbosity ( unsigned int  lvl )

inline

Set the verbosity of the output: 0 = silent, 1 = info, 2 = debug.

Definition at line 29 of file DTMeanTimerFitter.h.

References theVerbosityLevel.

{ theVerbosityLevel = lvl; }

Member Data Documentation

hDebugFile

TFile* DTMeanTimerFitter::hDebugFile

private

Definition at line 36 of file DTMeanTimerFitter.h.

Referenced by DTMeanTimerFitter(), evaluateVDriftAndReso(), and ~DTMeanTimerFitter().

hInputFile

TFile* DTMeanTimerFitter::hInputFile

private

Definition at line 37 of file DTMeanTimerFitter.h.

Referenced by evaluateVDriftAndReso().

theVerbosityLevel

unsigned int DTMeanTimerFitter::theVerbosityLevel

private

Definition at line 39 of file DTMeanTimerFitter.h.

Referenced by setVerbosity().