CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TH2Analyzer.cc
Go to the documentation of this file.
1 #include "Validation/RecoParticleFlow/interface/TH2Analyzer.h"
2 
3 #include "TH2D.h"
4 #include "TH1D.h"
5 #include "TF1.h"
6 #include "TPad.h"
7 
8 #include <string>
9 #include <iostream>
10 #include <sstream>
11 #include <algorithm>
12 
13 using namespace std;
14 
15 // remove?
16 void TH2Analyzer::Eval( int rebinFactor ) {
17  //std::cout << "Eval!" << std::endl;
18 
19  Reset();
20 
21  const string bname = hist2D_->GetName();
22  const string rebinName = bname + "_rebin";
23  rebinnedHist2D_ = (TH2D*) hist2D_ -> Clone( rebinName.c_str() );
24  rebinnedHist2D_->RebinX( rebinFactor );
25 
26  const string averageName = bname + "_average";
27  average_ = new TH1D( averageName.c_str(),"arithmetic average",
28  rebinnedHist2D_->GetNbinsX(),
29  rebinnedHist2D_->GetXaxis()->GetXmin(),
30  rebinnedHist2D_->GetXaxis()->GetXmax() );
31 
32  const string rmsName = bname + "_RMS";
33  RMS_ = new TH1D( rmsName.c_str(), "RMS",
34  rebinnedHist2D_->GetNbinsX(),
35  rebinnedHist2D_->GetXaxis()->GetXmin(),
36  rebinnedHist2D_->GetXaxis()->GetXmax() );
37 
38  const string sigmaGaussName = bname + "_sigmaGauss";
39  sigmaGauss_ = new TH1D(sigmaGaussName.c_str(), "sigmaGauss",
40  rebinnedHist2D_->GetNbinsX(),
41  rebinnedHist2D_->GetXaxis()->GetXmin(),
42  rebinnedHist2D_->GetXaxis()->GetXmax() );
43 
44  const string meanXName = bname + "_meanX";
45  meanXslice_ = new TH1D(meanXName.c_str(), "meanX",
46  rebinnedHist2D_->GetNbinsX(),
47  rebinnedHist2D_->GetXaxis()->GetXmin(),
48  rebinnedHist2D_->GetXaxis()->GetXmax() );
49 
50  ProcessSlices( rebinnedHist2D_ );
51 }
52 
53 void TH2Analyzer::Reset() {
54  if ( rebinnedHist2D_ ) delete rebinnedHist2D_;
55  if ( average_ ) delete average_;
56  if ( RMS_ ) delete RMS_;
57  if ( sigmaGauss_ ) delete sigmaGauss_;
58  if ( meanXslice_ ) delete meanXslice_;
59 
60  //for( unsigned i=0; i<parameters_.size(); ++i) {
61  // delete parameters_[i];
62  //}
63 
64  //parameters_.clear();
65 }
66 
67 void TH2Analyzer::Eval(const int rebinFactor, const int binxmin,
68  const int binxmax, const bool cst_binning)
69 {
70  Reset();
71  const string bname = hist2D_->GetName();
72  const string rebinName = bname + "_rebin";
73 
74  if (binxmax>hist2D_->GetNbinsX())
75  {
76  std::cout << "Error: TH2Analyzer.cc: binxmax>hist2D_->GetNbinsX()" << std::endl;
77  return;
78  }
79 
80  if (cst_binning)
81  {
82  //std::cout << "hist2D_->GetXaxis()->GetBinLowEdge(" << binxmin << ") = " << hist2D_->GetXaxis()->GetBinLowEdge(binxmin) << std::endl;
83  //std::cout << "hist2D_->GetXaxis()->GetBinUpEdge(" << binxmax << ") = " << hist2D_->GetXaxis()->GetBinUpEdge(binxmax) << std::endl;
84  //std::cout << "hist2D_->GetNbinsY() = " << hist2D_->GetNbinsY() << std::endl;
85  //std::cout << "hist2D_->GetYaxis()->GetXmin() = " << hist2D_->GetYaxis()->GetXmin() << std::endl;
86  //std::cout << "hist2D_->GetYaxis()->GetXmax() = " << hist2D_->GetYaxis()->GetXmax() << std::endl;
87  rebinnedHist2D_ = new TH2D(rebinName.c_str(),"rebinned histo",
88  binxmax-binxmin+1,
89  hist2D_->GetXaxis()->GetBinLowEdge(binxmin),
90  hist2D_->GetXaxis()->GetBinUpEdge(binxmax),
91  hist2D_->GetNbinsY(),
92  hist2D_->GetYaxis()->GetXmin(),
93  hist2D_->GetYaxis()->GetXmax() );
94  for (int binyc=1;binyc<hist2D_->GetNbinsY()+1;++binyc)
95  {
96  for (int binxc=binxmin;binxc<binxmax+1;++binxc)
97  {
98  //std::cout << "hist2D_->GetBinContent(" << binxc << "," << binyc << ") = " << hist2D_->GetBinContent(binxc,binyc) << std::endl;
99  //std::cout << "hist2D_->GetBinError(" << binxc << "," << binyc << ") = " << hist2D_->GetBinError(binxc,binyc) << std::endl;
100  //std::cout << "binxc-binxmin+1 = " << binxc-binxmin+1 << std::endl;
101  rebinnedHist2D_->SetBinContent(binxc-binxmin+1,binyc,hist2D_->GetBinContent(binxc,binyc));
102  rebinnedHist2D_->SetBinError(binxc-binxmin+1,binyc,hist2D_->GetBinError(binxc,binyc));
103  }
104  }
105  rebinnedHist2D_->RebinX( rebinFactor );
106 
107  const string averageName = bname + "_average";
108  average_ = new TH1D( averageName.c_str(),"arithmetic average",
109  rebinnedHist2D_->GetNbinsX(),
110  rebinnedHist2D_->GetXaxis()->GetXmin(),
111  rebinnedHist2D_->GetXaxis()->GetXmax() );
112 
113  const string rmsName = bname + "_RMS";
114  RMS_ = new TH1D( rmsName.c_str(), "RMS",
115  rebinnedHist2D_->GetNbinsX(),
116  rebinnedHist2D_->GetXaxis()->GetXmin(),
117  rebinnedHist2D_->GetXaxis()->GetXmax() );
118 
119  const string sigmaGaussName = bname + "_sigmaGauss";
120  sigmaGauss_ = new TH1D(sigmaGaussName.c_str(), "sigmaGauss",
121  rebinnedHist2D_->GetNbinsX(),
122  rebinnedHist2D_->GetXaxis()->GetXmin(),
123  rebinnedHist2D_->GetXaxis()->GetXmax() );
124 
125  const string meanXName = bname + "_meanX";
126  meanXslice_ = new TH1D(meanXName.c_str(), "meanX",
127  rebinnedHist2D_->GetNbinsX(),
128  rebinnedHist2D_->GetXaxis()->GetXmin(),
129  rebinnedHist2D_->GetXaxis()->GetXmax() );
130  }
131  else
132  {
133  // binning is not constant, but made to obtain almost the same number of events in each bin
134 
135  //std::cout << "binxmax-binxmin+1 = " << binxmax-binxmin+1 << std::endl;
136  //std::cout << "rebinFactor = " << rebinFactor << std::endl;
137  //std::cout << "(binxmax-binxmin+1)/rebinFactor = " << (binxmax-binxmin+1)/rebinFactor << std::endl;
138  //std::cout << "((binxmax-binxmin+1)%rebinFactor) = " << ((binxmax-binxmin+1)%rebinFactor) << std::endl;
139  //std::cout << "abs((binxmax-binxmin+1)/rebinFactor) = " << std::abs((binxmax-binxmin+1)/rebinFactor) << std::endl;
140 
141  unsigned int nbin=0;
142  if (((binxmax-binxmin+1)%rebinFactor)!=0.0)
143  {
144  nbin=std::abs((binxmax-binxmin+1)/rebinFactor)+1;
145  }
146  else nbin=(binxmax-binxmin+1)/rebinFactor;
147 
148  double *xlow = new double[nbin+1];
149  int *binlow = new int[nbin+1];
150 
151  TH1D* h0_slice1 = hist2D_->ProjectionY("h0_slice1", binxmin, binxmax, "");
152  const unsigned int totalNumberOfEvents=static_cast<unsigned int>(h0_slice1->GetEntries());
153  //std::cout << "totalNumberOfEvents = " << totalNumberOfEvents << std::endl;
154  delete h0_slice1;
155 
156  unsigned int neventsc=0;
157  unsigned int binXmaxc=binxmin+1;
158  xlow[0]=hist2D_->GetXaxis()->GetBinLowEdge(binxmin);
159  binlow[0]=binxmin;
160  for (unsigned int binc=1;binc<nbin;++binc)
161  {
162  while (neventsc<binc*totalNumberOfEvents/nbin)
163  {
164  TH1D* h0_slice1c = hist2D_->ProjectionY("h0_slice1",binxmin, binXmaxc, "");
165  neventsc=static_cast<unsigned int>(h0_slice1c->GetEntries());
166  // //std::cout << "FL : neventsc = " << neventsc << std::endl;
167  // //std::cout << "FL : binXmaxc = " << binXmaxc << std::endl;
168  ++binXmaxc;
169  delete h0_slice1c;
170  }
171  //std::cout << "binXmaxc-1 = " << binXmaxc-1 << std::endl;
172  binlow[binc]=binXmaxc-1;
173  xlow[binc]=hist2D_->GetXaxis()->GetBinLowEdge(binXmaxc-1);
174  }
175  xlow[nbin]=hist2D_->GetXaxis()->GetBinUpEdge(binxmax);
176  binlow[nbin]=binxmax;
177 
178  //for (unsigned int binc=0;binc<nbin+1;++binc)
179  //{
180  // std::cout << "xlow[" << binc << "] = " << xlow[binc] << std::endl;
181  //}
182 
183  rebinnedHist2D_ = new TH2D(rebinName.c_str(),"rebinned histo",
184  nbin, xlow, hist2D_->GetNbinsY(),
185  hist2D_->GetYaxis()->GetXmin(),
186  hist2D_->GetYaxis()->GetXmax() );
187  for (int binyc=1;binyc<hist2D_->GetNbinsY()+1;++binyc)
188  {
189  for (unsigned int binxc=1;binxc<nbin+1;++binxc)
190  {
191  double sum=0.0;
192  for (int binxh2c=binlow[binxc-1];binxh2c<binlow[binxc];++binxh2c)
193  {
194  sum+=hist2D_->GetBinContent(binxh2c,binyc);
195  }
196  rebinnedHist2D_->SetBinContent(binxc,binyc,sum);
197  }
198  }
199 
200  const string averageName = bname + "_average";
201  average_ = new TH1D( averageName.c_str(),"arithmetic average", nbin, xlow);
202 
203  const string rmsName = bname + "_RMS";
204  RMS_ = new TH1D( rmsName.c_str(), "RMS", nbin, xlow);
205 
206  const string sigmaGaussName = bname + "_sigmaGauss";
207  sigmaGauss_ = new TH1D(sigmaGaussName.c_str(), "sigmaGauss", nbin, xlow);
208 
209  const string meanXName = bname + "_meanX";
210  meanXslice_ = new TH1D(meanXName.c_str(), "meanX", nbin, xlow);
211  delete [] xlow;
212  delete [] binlow;
213  }
214  ProcessSlices( rebinnedHist2D_ );
215 }
216 
217 void TH2Analyzer::ProcessSlices( const TH2D* histo) {
218 
219  //std::cout << "ProcessSlices!" << std::endl;
220 
221  TH1::AddDirectory(0);
222 
223  for( int i=1; i<=histo->GetNbinsX(); ++i) {
224  TH1D* proj = histo->ProjectionY("toto", i, i);
225  const double mean = proj->GetMean();
226  const double rms = proj->GetRMS();
227  //std::cout << "mean = " << mean << std::endl;
228  //std::cout << "rms = " << rms << std::endl;
229  average_->SetBinContent( i, mean);
230  average_->SetBinError( i, proj->GetMeanError());
231  RMS_->SetBinContent(i, rms);
232  RMS_->SetBinError(i, proj->GetRMSError());
233 
234  const double error=(histo->GetXaxis()->GetBinUpEdge(i)-histo->GetXaxis()->GetBinLowEdge(i))/2.0;
235  //std::cout << "error = " << error << std::endl;
236  meanXslice_->SetBinContent(i, histo->GetXaxis()->GetBinLowEdge(i)+error);
237  meanXslice_->SetBinError(i, error);
238  //std::cout << "histo->GetXaxis()->GetBinLowEdge(" << i << ") = "
239  // << histo->GetXaxis()->GetBinLowEdge(i) << std::endl;
240  //std::cout << "meanXslice_->GetBinError(" << i << ") = "
241  // << meanXslice_->GetBinError(i) << std::endl;
242  ProcessSlice(i, proj );
243  delete proj;
244  }
245 
246  TH1::AddDirectory(1);
247 }
248 
249 
250 void TH2Analyzer::ProcessSlice(const int i, TH1D* proj ) const {
251 
252  //const double mean = proj->GetMean();
253  const double rms = proj->GetRMS();
254  //std::cout << "FL: mean = " << mean << std::endl;
255  //std::cout << "FL: rms = " << rms << std::endl;
256 
257  if (rms!=0.0)
258  {
259  const double fitmin=proj->GetMean()-proj->GetRMS();
260  const double fitmax=proj->GetMean()+proj->GetRMS();
261 
262  //std::cout << "i = " << i << std::endl;
263  //std::cout << "fitmin = " << fitmin << std::endl;
264  //std::cout << "fitmax = " << fitmax << std::endl;
265 
266  //proj->Draw();
267  TF1* f1= new TF1("f1", "gaus", fitmin, fitmax);
268  f1->SetParameters(proj->GetRMS(),proj->GetMean(),proj->GetBinContent(proj->GetMaximumBin()));
269  proj->Fit("f1","R", "", proj->GetXaxis()->GetXmin(), proj->GetXaxis()->GetXmax());
270 
271  //std::ostringstream oss;
272  //oss << i;
273  //const std::string plotfitname="Plots/fitbin_"+oss.str()+".eps";
274  //gPad->SaveAs( plotfitname.c_str() );
275  //std::cout << "param1 = " << f1->GetParameter(1) << std::endl;
276  //std::cout << "param2 = " << f1->GetParameter(2) << std::endl;
277  //std::cout << "paramError2 = " << f1->GetParError(2) << std::endl;
278 
279  sigmaGauss_->SetBinContent(i, f1->GetParameter(2));
280  sigmaGauss_->SetBinError(i, f1->GetParError(2));
281  delete f1;
282  }
283  else
284  {
285  sigmaGauss_->SetBinContent(i, 0.0);
286  sigmaGauss_->SetBinError(i, 0.0);
287  }
288 }
i
int i
Definition: DBlmapReader.cc:9
abs
#define abs(x)
Definition: mlp_lapack.h:159
TH2Analyzer::Eval
void Eval(const int rebinFactor)
Definition: TH2Analyzer.cc:16
TH2Analyzer::ProcessSlice
void ProcessSlice(const int i, TH1D *histo) const
Definition: TH2Analyzer.cc:250
TH2Analyzer::ProcessSlices
void ProcessSlices(const TH2D *histo)
Definition: TH2Analyzer.cc:217
TH2Analyzer::Reset
void Reset()
Definition: TH2Analyzer.cc:53
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
HcalObjRepresent::Reset
void Reset(std::vector< TH2F > &depth)
Definition: HcalObjRepresent.h:620
Generated for CMSSW Reference Manual by doxygen 1.8.5