CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DTMeanTimerFitter.cc
Go to the documentation of this file.
1 
2 /*
3  * See header file for a description of this class.
4  *
5  * $Date: 2009/02/26 09:46:57 $
6  * $Revision: 1.8 $
7  * \author S. Bolognesi - INFN Torino
8  */
9 
11 //#include "CalibMuon/DTCalibration/plugins/vDriftHistos.h"
13 
15 
16 #include <iostream>
17 
18 #include "TFile.h"
19 #include "TF1.h"
20 
21 using namespace std;
22 
23 DTMeanTimerFitter::DTMeanTimerFitter(TFile *file) : hInputFile(file), theVerbosityLevel(0) {
24  //hInputFile = new TFile("DTTMaxHistos.root", "READ");
25  hDebugFile = new TFile("DTMeanTimerFitter.root", "RECREATE");
26 }
27 
29  hDebugFile->Close();
30 }
31 
33 
34  // Retrieve histogram sets
36  vector<float> vDriftAndReso;
37 
38  // Check that the histo for this cell exists
39  if(histos->hTmax123 != 0) {
40  vector<TH1F*> hTMax; // histograms for <T_max> calculation
41  vector <TH1F*> hT0; // histograms for T0 evaluation
42  hTMax.push_back(histos->hTmax123);
43  hTMax.push_back(histos->hTmax124s72);
44  hTMax.push_back(histos->hTmax124s78);
45  hTMax.push_back(histos->hTmax134s72);
46  hTMax.push_back(histos->hTmax134s78);
47  hTMax.push_back(histos->hTmax234);
48 
49  hT0.push_back(histos->hTmax_3t0);
50  hT0.push_back(histos->hTmax_2t0);
51  hT0.push_back(histos->hTmax_t0);
52  hT0.push_back(histos->hTmax_0);
53 
54  vector<Double_t> factor; // factor relating the width of the Tmax distribution
55  // and the cell resolution
56  factor.push_back(sqrt(2./3.)); // hTmax123
57  factor.push_back(sqrt(2./7.)); // hTmax124s72
58  factor.push_back(sqrt(8./7.)); // hTmax124s78
59  factor.push_back(sqrt(2./7.)); // hTmax134s72
60  factor.push_back(sqrt(8./7.)); // hTmax134s78
61  factor.push_back(sqrt(2./3.)); // hTmax234
62 
63 
64  // Retrieve the gaussian mean and sigma for each TMax histogram
65  vector<Double_t> mean;
66  vector<Double_t> sigma;
67  vector<Double_t> count; //number of entries
68 
69  for(vector<TH1F*>::const_iterator ith = hTMax.begin();
70  ith != hTMax.end(); ith++) {
71  TF1 *funct = fitTMax(*ith);
72  if(!funct){
73  edm::LogError("DTMeanTimerFitter") << "Error when fitting TMax..histogram name" << (*ith)->GetName();
74  // return empty or -1 filled vector?
75  vector<float> defvec(6,-1);
76  return defvec;
77  }
78  hDebugFile->cd();
79  (*ith)->Write();
80 
81  // Get mean, sigma and number of entries of each histogram
82  mean.push_back(funct->GetParameter(1));
83  sigma.push_back(funct->GetParameter(2));
84  count.push_back((*ith)->GetEntries());
85  }
86 
87  Double_t tMaxMean=0.;
88  Double_t wTMaxSum=0.;
89  Double_t sigmaT=0.;
90  Double_t wSigmaSum = 0.;
91 
92  //calculate total mean and sigma
93  for(int i=0; i<=5; i++) {
94  if(count[i]<200) continue;
95  tMaxMean += mean[i]*(count[i]/(sigma[i]*sigma[i]));
96  wTMaxSum += count[i]/(sigma[i]*sigma[i]);
97  sigmaT += count[i]*factor[i]*sigma[i];
98  wSigmaSum += count[i];
99  // cout << "TMaxMean "<<i<<": "<< mean[i] << " entries: " << count[i]
100  // << " sigma: " << sigma[i]
101  // << " weight: " << (count[i]/(sigma[i]*sigma[i])) << endl;
102  }
103  if((!wTMaxSum)||(!wSigmaSum)){
104  edm::LogError("DTMeanTimerFitter") << "Error zero sum of weights..returning default";
105  vector<float> defvec(6,-1);
106  return defvec;
107  }
108 
109  tMaxMean /= wTMaxSum;
110  sigmaT /= wSigmaSum;
111 
112  //calculate v_drift and resolution
113  Double_t vDrift = 2.1 / tMaxMean; //2.1 is the half cell length in cm
114  Double_t reso = vDrift * sigmaT;
115  vDriftAndReso.push_back(vDrift);
116  vDriftAndReso.push_back(reso);
117  //if(theVerbosityLevel >= 1)
118  edm::LogVerbatim("DTMeanTimerFitter") << " final TMaxMean=" << tMaxMean << " sigma= " << sigmaT
119  << " v_d and reso: " << vDrift << " " << reso << endl;
120 
121  // Order t0 histogram by number of entries (choose histograms with higher nr. of entries)
122  map<Double_t,TH1F*> hEntries;
123  for(vector<TH1F*>::const_iterator ith = hT0.begin();
124  ith != hT0.end(); ith++) {
125  hEntries[(*ith)->GetEntries()] = (*ith);
126  }
127 
128  // add at the end of hT0 the two hists with the higher number of entries
129  int counter = 0;
130  for(map<Double_t,TH1F*>::reverse_iterator iter = hEntries.rbegin();
131  iter != hEntries.rend(); iter++) {
132  counter++;
133  if (counter==1) hT0.push_back(iter->second);
134  else if (counter==2) {hT0.push_back(iter->second); break;}
135  }
136 
137  // Retrieve the gaussian mean and sigma of histograms for Delta(t0) evaluation
138  vector<Double_t> meanT0;
139  vector<Double_t> sigmaT0;
140  vector<Double_t> countT0;
141 
142  for(vector<TH1F*>::const_iterator ith = hT0.begin();
143  ith != hT0.end(); ith++) {
144  try{
145  (*ith)->Fit("gaus");
146  } catch(...){
147  edm::LogError("DTMeanTimerFitter") << "Exception when fitting T0..histogram " << (*ith)->GetName();
148  // return empty or -1 filled vector?
149  vector<float> defvec(6,-1);
150  return defvec;
151  }
152  TF1 *f1 = (*ith)->GetFunction("gaus");
153  // Get mean, sigma and number of entries of the histograms
154  meanT0.push_back(f1->GetParameter(1));
155  sigmaT0.push_back(f1->GetParameter(2));
156  countT0.push_back((*ith)->GetEntries());
157  }
158  //calculate Delta(t0)
159  if(hT0.size() != 6) { // check if you have all the t0 hists
160  edm::LogVerbatim("DTMeanTimerFitter") << "t0 histograms = " << hT0.size();
161  for(int i=1; i<=4;i++) {
162  vDriftAndReso.push_back(-1);
163  }
164  return vDriftAndReso;
165  }
166 
167  for(int it0=0; it0<=2; it0++) {
168  if((countT0[it0] > 200) && (countT0[it0+1] > 200)) {
169  Double_t deltaT0 = meanT0[it0] - meanT0[it0+1];
170  vDriftAndReso.push_back(deltaT0);
171  }
172  else
173  vDriftAndReso.push_back(999.);
174  }
175  //deltat0 using hists with max nr. of entries
176  if((countT0[4] > 200) && (countT0[5] > 200)) {
177  Double_t t0Diff = histos->GetT0Factor(hT0[4]) - histos->GetT0Factor(hT0[5]);
178  Double_t deltaT0MaxEntries = (meanT0[4] - meanT0[5])/ t0Diff;
179  vDriftAndReso.push_back(deltaT0MaxEntries);
180  }
181  else
182  vDriftAndReso.push_back(999.);
183  }
184  else {
185  for(int i=1; i<=6; i++) {
186  // 0=vdrift, 1=reso, 2=(3deltat0-2deltat0), 3=(2deltat0-1deltat0),
187  // 4=(1deltat0-0deltat0), 5=deltat0 from hists with max entries,
188  vDriftAndReso.push_back(-1);
189  }
190  }
191  return vDriftAndReso;
192 }
193 
194 
195 
197  // Find distribution peak and fit range
198  Double_t peak = (((((histo->GetXaxis())->GetXmax())-((histo->GetXaxis())->GetXmin()))/histo->GetNbinsX())*
199  (histo->GetMaximumBin()))+((histo->GetXaxis())->GetXmin());
200  //if(theVerbosityLevel >= 1)
201  LogDebug("DTMeanTimerFitter") <<"Peak "<<peak<<" : "<<"xmax "<<((histo->GetXaxis())->GetXmax())
202  <<" xmin "<<((histo->GetXaxis())->GetXmin())
203  <<" nbin "<<histo->GetNbinsX()
204  <<" bin with max "<<(histo->GetMaximumBin());
205  Double_t range = 2.*histo->GetRMS();
206 
207  // Fit each Tmax histogram with a Gaussian in a restricted interval
208  TF1 *rGaus = new TF1("rGaus","gaus",peak-range,peak+range);
209  rGaus->SetMarkerSize();// just silence gcc complainining about unused var
210  try{
211  histo->Fit("rGaus","R");
212  } catch(...){
213  edm::LogError("DTMeanTimerFitter") << "Exception when fitting TMax..histogram " << histo->GetName()
214  << " setting return function pointer to zero";
215  return 0;
216  }
217  TF1 *f1 = histo->GetFunction("rGaus");
218  return f1;
219  }
#define LogDebug(id)
int i
Definition: DBlmapReader.cc:9
TH1F * hTmax134s78
Definition: vDriftHistos.h:410
TH1F * hTmax_2t0
Definition: vDriftHistos.h:419
TH1F * hTmax_3t0
Definition: vDriftHistos.h:412
virtual ~DTMeanTimerFitter()
Destructor.
DTMeanTimerFitter(TFile *file)
Constructor.
int GetT0Factor(TH1F *hist)
Definition: vDriftHistos.h:391
TH1F * hTmax234
Definition: vDriftHistos.h:411
TH1F * hTmax_t0
Definition: vDriftHistos.h:426
TH1F * hTmax134s72
Definition: vDriftHistos.h:409
TH1F * hTmax124s72
Definition: vDriftHistos.h:407
TF1 * fitTMax(TH1F *histo)
Really do the fit.
TH1F * hTmax_0
Definition: vDriftHistos.h:433
std::vector< float > evaluateVDriftAndReso(TString N)
Fit the TMax histos and evaluate VDrift and resolution.
T sqrt(T t)
Definition: SSEVec.h:46
TH1F * hTmax124s78
Definition: vDriftHistos.h:408
#define N
Definition: blowfish.cc:9
TH1F * hTmax123
Definition: vDriftHistos.h:406