DTCalibration/src/DTMeanTimerFitter.cc

0001
0002 /*
0003  *  See header file for a description of this class.
0004  *
0005  *  \author S. Bolognesi - INFN Torino
0006  */
0007
0008 #include "CalibMuon/DTCalibration/interface/DTMeanTimerFitter.h"
0009 //#include "CalibMuon/DTCalibration/plugins/vDriftHistos.h"
0010 #include "CalibMuon/DTCalibration/interface/vDriftHistos.h"
0011
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include <iostream>
0015
0016 #include "TFile.h"
0017 #include "TF1.h"
0018
0019 using namespace std;
0020
0021 DTMeanTimerFitter::DTMeanTimerFitter(TFile *file) : hInputFile(file), theVerbosityLevel(0) {
0022   //hInputFile = new TFile("DTTMaxHistos.root", "READ");
0023   hDebugFile = new TFile("DTMeanTimerFitter.root", "RECREATE");
0024 }
0025
0026 DTMeanTimerFitter::~DTMeanTimerFitter() { hDebugFile->Close(); }
0027
0028 vector<float> DTMeanTimerFitter::evaluateVDriftAndReso(const TString &N) {
0029   // Retrieve histogram sets
0030   hTMaxCell *histos = new hTMaxCell(N, hInputFile);
0031   vector<float> vDriftAndReso;
0032
0033   // Check that the histo for this cell exists
0034   if (histos->hTmax123 != nullptr) {
0035     vector<TH1F *> hTMax;  // histograms for <T_max> calculation
0036     vector<TH1F *> hT0;    // histograms for T0 evaluation
0037     hTMax.push_back(histos->hTmax123);
0038     hTMax.push_back(histos->hTmax124s72);
0039     hTMax.push_back(histos->hTmax124s78);
0040     hTMax.push_back(histos->hTmax134s72);
0041     hTMax.push_back(histos->hTmax134s78);
0042     hTMax.push_back(histos->hTmax234);
0043
0044     hT0.push_back(histos->hTmax_3t0);
0045     hT0.push_back(histos->hTmax_2t0);
0046     hT0.push_back(histos->hTmax_t0);
0047     hT0.push_back(histos->hTmax_0);
0048
0049     vector<Double_t> factor;          // factor relating the width of the Tmax distribution
0050                                       // and the cell resolution
0051     factor.push_back(sqrt(2. / 3.));  // hTmax123
0052     factor.push_back(sqrt(2. / 7.));  // hTmax124s72
0053     factor.push_back(sqrt(8. / 7.));  // hTmax124s78
0054     factor.push_back(sqrt(2. / 7.));  // hTmax134s72
0055     factor.push_back(sqrt(8. / 7.));  // hTmax134s78
0056     factor.push_back(sqrt(2. / 3.));  // hTmax234
0057
0058     // Retrieve the gaussian mean and sigma for each TMax histogram
0059     vector<Double_t> mean;
0060     vector<Double_t> sigma;
0061     vector<Double_t> count;  //number of entries
0062
0063     for (vector<TH1F *>::const_iterator ith = hTMax.begin(); ith != hTMax.end(); ++ith) {
0064       TF1 *funct = fitTMax(*ith);
0065       if (!funct) {
0066         edm::LogError("DTMeanTimerFitter") << "Error when fitting TMax..histogram name" << (*ith)->GetName();
0067         // return empty or -1 filled vector?
0068         vector<float> defvec(6, -1);
0069         return defvec;
0070       }
0071       hDebugFile->cd();
0072       (*ith)->Write();
0073
0074       // Get mean, sigma and number of entries of each histogram
0075       mean.push_back(funct->GetParameter(1));
0076       sigma.push_back(funct->GetParameter(2));
0077       count.push_back((*ith)->GetEntries());
0078     }
0079
0080     Double_t tMaxMean = 0.;
0081     Double_t wTMaxSum = 0.;
0082     Double_t sigmaT = 0.;
0083     Double_t wSigmaSum = 0.;
0084
0085     //calculate total mean and sigma
0086     for (int i = 0; i <= 5; i++) {
0087       if (count[i] < 200)
0088         continue;
0089       tMaxMean += mean[i] * (count[i] / (sigma[i] * sigma[i]));
0090       wTMaxSum += count[i] / (sigma[i] * sigma[i]);
0091       sigmaT += count[i] * factor[i] * sigma[i];
0092       wSigmaSum += count[i];
0093       // cout << "TMaxMean "<<i<<": "<< mean[i] << " entries: " << count[i]
0094       // << " sigma: " << sigma[i]
0095       // << " weight: " << (count[i]/(sigma[i]*sigma[i])) << endl;
0096     }
0097     if ((!wTMaxSum) || (!wSigmaSum)) {
0098       edm::LogError("DTMeanTimerFitter") << "Error zero sum of weights..returning default";
0099       vector<float> defvec(6, -1);
0100       return defvec;
0101     }
0102
0103     tMaxMean /= wTMaxSum;
0104     sigmaT /= wSigmaSum;
0105
0106     //calculate v_drift and resolution
0107     Double_t vDrift = 2.1 / tMaxMean;  //2.1 is the half cell length in cm
0108     Double_t reso = vDrift * sigmaT;
0109     vDriftAndReso.push_back(vDrift);
0110     vDriftAndReso.push_back(reso);
0111     //if(theVerbosityLevel >= 1)
0112     edm::LogVerbatim("DTMeanTimerFitter")
0113         << " final TMaxMean=" << tMaxMean << " sigma= " << sigmaT << " v_d and reso: " << vDrift << " " << reso << endl;
0114
0115     // Order t0 histogram by number of entries (choose histograms with higher nr. of entries)
0116     map<Double_t, TH1F *> hEntries;
0117     for (vector<TH1F *>::const_iterator ith = hT0.begin(); ith != hT0.end(); ++ith) {
0118       hEntries[(*ith)->GetEntries()] = (*ith);
0119     }
0120
0121     // add at the end of hT0 the two hists with the higher number of entries
0122     int counter = 0;
0123     for (map<Double_t, TH1F *>::reverse_iterator iter = hEntries.rbegin(); iter != hEntries.rend(); ++iter) {
0124       counter++;
0125       if (counter == 1)
0126         hT0.push_back(iter->second);
0127       else if (counter == 2) {
0128         hT0.push_back(iter->second);
0129         break;
0130       }
0131     }
0132
0133     // Retrieve the gaussian mean and sigma of histograms for Delta(t0) evaluation
0134     vector<Double_t> meanT0;
0135     vector<Double_t> sigmaT0;
0136     vector<Double_t> countT0;
0137
0138     for (vector<TH1F *>::const_iterator ith = hT0.begin(); ith != hT0.end(); ++ith) {
0139       try {
0140         (*ith)->Fit("gaus");
0141       } catch (std::exception const &) {
0142         edm::LogError("DTMeanTimerFitter") << "Exception when fitting T0..histogram " << (*ith)->GetName();
0143         // return empty or -1 filled vector?
0144         vector<float> defvec(6, -1);
0145         return defvec;
0146       }
0147       TF1 *f1 = (*ith)->GetFunction("gaus");
0148       // Get mean, sigma and number of entries of the  histograms
0149       meanT0.push_back(f1->GetParameter(1));
0150       sigmaT0.push_back(f1->GetParameter(2));
0151       countT0.push_back((*ith)->GetEntries());
0152     }
0153     //calculate Delta(t0)
0154     if (hT0.size() != 6) {  // check if you have all the t0 hists
0155       edm::LogVerbatim("DTMeanTimerFitter") << "t0 histograms = " << hT0.size();
0156       for (int i = 1; i <= 4; i++) {
0157         vDriftAndReso.push_back(-1);
0158       }
0159       return vDriftAndReso;
0160     }
0161
0162     for (int it0 = 0; it0 <= 2; it0++) {
0163       if ((countT0[it0] > 200) && (countT0[it0 + 1] > 200)) {
0164         Double_t deltaT0 = meanT0[it0] - meanT0[it0 + 1];
0165         vDriftAndReso.push_back(deltaT0);
0166       } else
0167         vDriftAndReso.push_back(999.);
0168     }
0169     //deltat0 using hists with max nr. of entries
0170     if ((countT0[4] > 200) && (countT0[5] > 200)) {
0171       Double_t t0Diff = histos->GetT0Factor(hT0[4]) - histos->GetT0Factor(hT0[5]);
0172       Double_t deltaT0MaxEntries = (meanT0[4] - meanT0[5]) / t0Diff;
0173       vDriftAndReso.push_back(deltaT0MaxEntries);
0174     } else
0175       vDriftAndReso.push_back(999.);
0176   } else {
0177     for (int i = 1; i <= 6; i++) {
0178       // 0=vdrift, 1=reso,  2=(3deltat0-2deltat0), 3=(2deltat0-1deltat0),
0179       // 4=(1deltat0-0deltat0), 5=deltat0 from hists with max entries,
0180       vDriftAndReso.push_back(-1);
0181     }
0182   }
0183   return vDriftAndReso;
0184 }
0185
0186 TF1 *DTMeanTimerFitter::fitTMax(TH1F *histo) {
0187   // Find distribution peak and fit range
0188   Double_t peak = (((((histo->GetXaxis())->GetXmax()) - ((histo->GetXaxis())->GetXmin())) / histo->GetNbinsX()) *
0189                    (histo->GetMaximumBin())) +
0190                   ((histo->GetXaxis())->GetXmin());
0191   //if(theVerbosityLevel >= 1)
0192   LogDebug("DTMeanTimerFitter") << "Peak " << peak << " : "
0193                                 << "xmax " << ((histo->GetXaxis())->GetXmax()) << "            xmin "
0194                                 << ((histo->GetXaxis())->GetXmin()) << "            nbin " << histo->GetNbinsX()
0195                                 << "            bin with max " << (histo->GetMaximumBin());
0196   Double_t range = 2. * histo->GetRMS();
0197
0198   // Fit each Tmax histogram with a Gaussian in a restricted interval
0199   TF1 *rGaus = new TF1("rGaus", "gaus", peak - range, peak + range);
0200   rGaus->SetMarkerSize();  // just silence gcc complainining about unused var
0201   try {
0202     histo->Fit("rGaus", "R");
0203   } catch (std::exception const &) {
0204     edm::LogError("DTMeanTimerFitter") << "Exception when fitting TMax..histogram " << histo->GetName()
0205                                        << "   setting return function pointer to zero";
0206     return nullptr;
0207   }
0208   TF1 *f1 = histo->GetFunction("rGaus");
0209   return f1;
0210 }