CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
fitMaximumPoint.cc File Reference

Go to the source code of this file.

Functions

Float_t Fit_MaximumPoint (TH2F *h2, Int_t Ireturn=0)
 
Double_t Pol2 (Double_t *x, Double_t *par)
 
Double_t Pol2_Special (Double_t *x, Double_t *par)
 

Function Documentation

Float_t Fit_MaximumPoint ( TH2F *  h2,
Int_t  Ireturn = 0 
)

Definition at line 24 of file fitMaximumPoint.cc.

References a, b, EnergyCorrector::c, gather_cfg::cout, pyrootRender::da, EcalCondDB::db, python.connectstrParser::f1, fitWZ::gMinuit, i, Pol2(), Pol2_Special(), mathSSE::sqrt(), and MetTreeProducer::var().

24  { // xxx
25 //====================================================
26 //====================================================
27 
28  // -------------------------------------------
29  // 1) Clean beam profile
30  // o all bins: only bins > 2 entries
31  // o per x-bin: only bins > 75% of maximum
32  // -------------------------------------------
33  // get some caracteristics
34  Int_t x_nbins = h2->GetXaxis()->GetNbins();
35  Float_t x_max = h2->GetXaxis()->GetXmax();
36  Float_t x_min = h2->GetXaxis()->GetXmin();
37  Float_t delta_x = h2->GetXaxis()->GetBinWidth(10);
38  Int_t y_nbins = h2->GetYaxis()->GetNbins();
39  Float_t y_max = h2->GetYaxis()->GetXmax();
40  Float_t y_min = h2->GetYaxis()->GetXmin();
41  Float_t delta_y = h2->GetYaxis()->GetBinWidth(10);
42  TH2F *IHBeamProf_clean_1 = new TH2F("IHBeamProf_clean_1", "Max sam versus X (clean) ", x_nbins, x_min, x_max, y_nbins, y_min, y_max);
43  TH2F *IHBeamProf_clean_2 = new TH2F("IHBeamProf_clean_2", "Max sam versus X (clean) ", x_nbins, x_min, x_max, y_nbins, y_min, y_max);
44 
45  // [1a] Only keep bins with more than 2 entries
46  // Also remember for each x-bin the maximum bin in y that satisfies this result
47  Int_t y_binmax_high[500] = {1};
48  for (Int_t ix=1;ix<x_nbins+1;ix++){
49  for (Int_t iy=1;iy<y_nbins+1;iy++){
50  Int_t Nevts_clean = 0;
51  if( h2->GetBinContent(ix,iy) > 1 ) {
52  IHBeamProf_clean_1->SetBinContent(ix,iy, h2->GetBinContent(ix,iy));
53  if( iy > y_binmax_high[ix] ){ y_binmax_high[ix] = iy; }
54  }
55  }
56  }
57  // [1b] Only keep events with more than 85% of the maximum
58  Int_t y_binmax_low[500] = {1};
59  for (Int_t ix=1;ix<x_nbins+1;ix++){
60  y_binmax_low[ix] = (Int_t)( 0.85*(Float_t)( y_binmax_high[ix]));
61  for (Int_t iy=y_binmax_low[ix];iy<y_binmax_high[ix]+1;iy++){
62  IHBeamProf_clean_2->SetBinContent(ix,iy,IHBeamProf_clean_1->GetBinContent(ix,iy));
63  }
64  }
65 
66 
67  // -----------------------------------------------------------------------
68  // 2) Find region to fit
69  // o Make profile by compuing means in every bin (error = RMS/sqrt(N))
70  // o Find maximum
71  // -----------------------------------------------------------------------
72  TH1F *IHBeamProf_clean_3 = new TH1F((h2->GetName()+(TString("_clprof"))).Data(), "Meam Energy versus X (clean) ", x_nbins, x_min, x_max);
73 
74  // [2a] Make TH1F containing the mean values
75  TH1F *h_temp = new TH1F("h_temp"," Energy distribution for single bin in X", y_nbins, y_min, y_max);
76  for (Int_t ix=1;ix<x_nbins+1;ix++){
77  double Nevt_slice = 0;
78  //Int_t Nevt_slice = 0;
79  for (Int_t iy=1;iy<y_nbins+1;iy++){
80  Int_t Nevt_tmp = (int)IHBeamProf_clean_2->GetBinContent(ix,iy);
81  Nevt_slice += Nevt_tmp;
82  h_temp->SetBinContent(iy,Nevt_tmp);
83  }
84  Float_t Y_mean = h_temp->GetMean();
85  Float_t Y_rms = h_temp->GetRMS();
86  Float_t Y_mean_error = (Nevt_slice>0) ? Y_rms/sqrt(Nevt_slice) : 9999.;
87  IHBeamProf_clean_3->SetBinContent(ix,Y_mean);
88  IHBeamProf_clean_3->SetBinError(ix,Y_mean_error);
89  printf("%d %f %f %d %f\n",ix,Y_mean,Y_rms,Nevt_slice,Y_mean_error);
90  }
91 
92  // [2b] Find maximum
93  Float_t Y_tmp_max = -999.;
94  Float_t x_max_guess = 0.;
95  Int_t x_max_guess_bin = 0;
96  for (Int_t ix=1;ix<x_nbins+1;ix++){
97  Float_t X_tmp = IHBeamProf_clean_3->GetBinCenter(ix);
98  Float_t Y_tmp = IHBeamProf_clean_3->GetBinContent(ix);
99  Float_t Y_tmp_error = IHBeamProf_clean_3->GetBinError(ix);
100  if( (fabs(X_tmp) < 10.) && (Y_tmp > Y_tmp_max) && (Y_tmp_error < 100.)){
101  Y_tmp_max = Y_tmp;
102  x_max_guess = X_tmp;
103  x_max_guess_bin = ix;
104  printf("Xtmp = %5.3f Ytmp = %5.3f Ytmp_max = %5.3f\n",X_tmp,Y_tmp,Y_tmp_max);
105  }
106  }
107  printf("Fit: o Start for maximum = %5.3f\n",x_max_guess);
108 
109  // [2c] Define the fit range (0.975% and error should be less than 8\%)
110  Int_t fitbinmin = 9999;
111  Int_t fitbinmax = -9999;
112  Float_t x_tmp,y_tmp,e_tmp,er_tmp;
113  for (Int_t ix = x_max_guess_bin; ix<x_nbins+1;ix++){
114  fitbinmax = ix;
115  x_tmp = IHBeamProf_clean_3->GetBinCenter(ix);
116  y_tmp = IHBeamProf_clean_3->GetBinContent(ix);
117  e_tmp = IHBeamProf_clean_3->GetBinError(ix);
118  er_tmp = (y_tmp>0.) ? (e_tmp/y_tmp) : 1.00;
119  printf("%3d %f %f %f -- %f %f\n",ix,x_tmp,y_tmp,e_tmp,y_tmp/Y_tmp_max,er_tmp);
120  if( y_tmp < 0.975*Y_tmp_max && er_tmp < 0.008) break;
121  }
122  for (Int_t ix = x_max_guess_bin; ix>1 ;ix--){
123  fitbinmin = ix;
124  x_tmp = IHBeamProf_clean_3->GetBinCenter(ix);
125  y_tmp = IHBeamProf_clean_3->GetBinContent(ix);
126  e_tmp = IHBeamProf_clean_3->GetBinError(ix);
127  er_tmp = (y_tmp>0.) ? (e_tmp/y_tmp) : 1.00;
128  printf("%3d %f %f %f -- %f %f\n",ix,x_tmp,y_tmp,e_tmp,y_tmp/Y_tmp_max,er_tmp);
129  if( y_tmp < 0.975* Y_tmp_max && er_tmp < 0.008) break;
130  }
131 
132  Double_t posmin_fit = x_min+fitbinmin*delta_x;
133  Double_t posmax_fit = x_min+fitbinmax*delta_x;
134  printf(" o Fit range = %5.3f -- %5.3f -- %5.3f\n",posmin_fit,x_max_guess,posmax_fit);
135  if( fabs(posmax_fit - posmin_fit ) < 4. ){
136  printf("Something is wrong with this range: returning dummy value\n");
137  posmin_fit = x_max_guess - 6.0;
138  posmax_fit = x_max_guess + 6.0;
139  return -99.00;
140  }
141 
142  // -------------------------------------------------
143  // 3) Do the actual fit
144  // o make clone of histogram
145  // -------------------------------------------------
146  TH1F *h1 = (TH1F *) IHBeamProf_clean_3->Clone();
147  // 3a] Do the actual fit
148  Double_t fitresults[3] = {0.};
149  Double_t fitresults_errors[3] = {0.};
150  TF1 *f1 = new TF1("f1",Pol2,-50.,50.,3);
151  f1->SetParameters( 1.,1.,1.);
152  h1->Fit(f1,"Q","",posmin_fit, posmax_fit);
153  for(int i=0 ; i< 3 ; i++) {
154  fitresults[i] = f1->GetParameter(i);
155  fitresults_errors[i] = f1->GetParError(i);
156  }
157  Float_t chi2 = f1->GetChisquare()/f1->GetNDF();
158  Float_t a = fitresults[2]; Float_t da = fitresults_errors[2];
159  Float_t b = fitresults[1]; Float_t db = fitresults_errors[1];
160  Float_t c = fitresults[0]; Float_t dc = fitresults_errors[0];
161  Float_t x0 = (-1.*b)/(2*a);
162  printf(" a = %7.2f b = %7.2f c = %7.2f\n",a,b,c);
163  printf(" da = %7.2f db = %7.2f dc = %7.2f\n",da,db,dc);
164 
165  cout << "risultati del fit polinomiale: " << fitresults[0] << " " << fitresults[1] << " " << fitresults[2] << endl;
166 
167  char myProfTitle[200];
168  sprintf(myProfTitle, h2->GetTitle());
169  strcat(myProfTitle,"_prof");
170  h1->Write(myProfTitle);
171  //h1->Write(); //write the profile
172  // ----------------------------
173  // 4) Compute uncertainty on x0
174  // ----------------------------
175  // [4a] compute dxo using the covariance matrix
176  // covariance matrix
177  Double_t CovMat[3][3];
178  gMinuit->mnemat(&CovMat[0][0],3);
179  Float_t v11 = CovMat[0][0]; Float_t v12 = CovMat[0][1]; Float_t v13 = CovMat[0][2];
180  Float_t v21 = CovMat[1][0]; Float_t v22 = CovMat[1][1]; Float_t v23 = CovMat[1][2];
181  Float_t v31 = CovMat[2][0]; Float_t v32 = CovMat[2][1]; Float_t v33 = CovMat[2][2];
182  printf("Covariance Matrix: v11 = %f v12 = %f v13 = %f\n",CovMat[0][0],CovMat[0][1],CovMat[0][2]);
183  printf(" v21 = %f v22 = %f v23 = %f\n",CovMat[1][0],CovMat[1][1],CovMat[1][2]);
184  printf(" v31 = %f v32 = %f v33 = %f\n",CovMat[2][0],CovMat[2][1],CovMat[2][2]);
185  // jacobiaan
186  Float_t j1 = b/(2*a*a);
187  Float_t j2 = -1./(2*a);
188  // determinant covariance martix
189  Float_t det = v11*(v33*v22-v32*v23)-v21*(v33*v12-v32*v13)+v31*(v23*v12-v22*v13);
190  printf("Determinant = %f\n",det);
191  // inverse matrix
192  Float_t v11_i = v33*v22-v32*v23; Float_t v12_i = -(v33*v12-v32*v13); Float_t v13_i = v23*v12-v22*v13;
193  Float_t v21_i = -(v33*v21-v31*v23); Float_t v22_i = v33*v11-v31*v13 ; Float_t v23_i = -(v23*v11-v21*v13) ;
194  Float_t v31_i = v32*v21-v31*v22; Float_t v32_i = -(v32*v11-v31*v12); Float_t v33_i = v22*v11-v21*v12;
195  // variance
196  Float_t var = j1*(j1*v11_i+j2*v12_i)+j2*(j1*v21_i+j2*v22_i);
197  var /= det;
198  Float_t dx0 = sqrt(var);
199  printf("Type 1 fit: o x0 = %f +/- %f\n",x0,dx0);
200 
201 
202  // ---------------------
203  // 5) Second type of fit
204  // ---------------------
205  // 5a] Do the actual fit
206  TH1F *h0 = (TH1F *) IHBeamProf_clean_3->Clone();
207  Double_t fitresults0[3] = {0.};
208  Double_t fitresults0_errors[3] = {0.};
209  TF1 *f0 = new TF1("f0",Pol2_Special,-50.,50.,3);
210  f0->SetParameters( -1.,x_max_guess,5000.);
211  h0->Fit(f0,"Q","",posmin_fit, posmax_fit);
212  for(int i=0 ; i< 3 ; i++) {
213  fitresults0[i] = f0->GetParameter(i);
214  fitresults0_errors[i] = f0->GetParError(i);
215  }
216  Float_t a0 = fitresults0[2]; Float_t da0 = fitresults0_errors[2];
217  Float_t b0 = fitresults0[1]; Float_t db0 = fitresults0_errors[1];
218  Float_t c0 = fitresults0[0]; Float_t dc0 = fitresults0_errors[0];
219 
220  Double_t CovMat0[3][3];
221  gMinuit->mnemat(&CovMat0[0][0],3);
222  //printf("Cov0[1][1] = %f\n",CovMat0[1][1]);
223  Float_t db0cov = sqrt(CovMat0[1][1]);
224  printf("Type 2 fit: o x0 = %f +/- %f %f \n",b0,db0,db0cov);
225  delete IHBeamProf_clean_1;
226  delete IHBeamProf_clean_2;
227  delete IHBeamProf_clean_3;
228  delete h_temp;
229  delete f1;
230  delete f0;
231 
232  if(Ireturn == 10) {return x0;}
233  if(Ireturn == 11) {return dx0;}
234  if(Ireturn == 20) {return b0;}
235  if(Ireturn == 21) {return db0cov;}
236  if(Ireturn == 30) {return chi2;}
237 
238  if(Ireturn == 99) {return c0;}
239 
240 
241 
242 }
i
int i
Definition: DBlmapReader.cc:9
Pol2
Double_t Pol2(Double_t *x, Double_t *par)
Definition: fitMaximumPoint.cc:1
EcalCondDB.db
tuple db
Definition: EcalCondDB.py:151
Pol2_Special
Double_t Pol2_Special(Double_t *x, Double_t *par)
Definition: fitMaximumPoint.cc:11
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
fitWZ.gMinuit
tuple gMinuit
Definition: fitWZ.py:35
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
Double_t Pol2 ( Double_t *  x,
Double_t *  par 
)

Definition at line 1 of file fitMaximumPoint.cc.

Referenced by Fit_MaximumPoint().

2 {
3  //polynoom(2)
4  Double_t part0 = par[0];
5  Double_t part1 = par[1]*x[0];
6  Double_t part2 = par[2]*x[0]*x[0];
7  Double_t fitval = part0 + part1 + part2;
8  return fitval;
9 }
x
T x() const
Cartesian x coordinate.
Definition: Basic3DVectorLD.h:127
Double_t Pol2_Special ( Double_t *  x,
Double_t *  par 
)

Definition at line 11 of file fitMaximumPoint.cc.

Referenced by Fit_MaximumPoint().

12 {
13  //polynoom(2)
14  Double_t part0 = par[0];
15  Double_t part1 = -2.*par[2]*par[1]*x[0];
16  Double_t part2 = par[2]*x[0]*x[0];
17  Double_t fitval = part0 + part1 + part2;
18  return fitval;
19 }
x
T x() const
Cartesian x coordinate.
Definition: Basic3DVectorLD.h:127
Generated for CMSSW Reference Manual by doxygen 1.8.5