CMS 3D CMS Logo
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

◆ Fit_MaximumPoint()

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

Definition at line 24 of file fitMaximumPoint.cc.

References a, a0, b, b0, c, fftjetpileupestimator_calo_uncalib_cfi::c0, hltPixelTracks_cff::chi2, gather_cfg::cout, Calorimetry_cff::da, dqmiodatasetharvest::db, DeadROC_duringRun::f1, mps_fire::i, createfilelist::int, Pol2(), Pol2_Special(), mathSSE::sqrt(), and trigObjTnPSource_cfi::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 }
Calorimetry_cff.da
da
Definition: Calorimetry_cff.py:155
mps_fire.i
i
Definition: mps_fire.py:428
hltPixelTracks_cff.chi2
chi2
Definition: hltPixelTracks_cff.py:25
trigObjTnPSource_cfi.var
var
Definition: trigObjTnPSource_cfi.py:21
Pol2
Double_t Pol2(Double_t *x, Double_t *par)
Definition: fitMaximumPoint.cc:1
DeadROC_duringRun.f1
f1
Definition: DeadROC_duringRun.py:219
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:19
a0
static constexpr float a0
Definition: L1EGammaCrystalsEmulatorProducer.cc:79
b
double b
Definition: hdecay.h:118
a
double a
Definition: hdecay.h:119
fftjetpileupestimator_calo_uncalib_cfi.c0
c0
Definition: fftjetpileupestimator_calo_uncalib_cfi.py:8
gather_cfg.cout
cout
Definition: gather_cfg.py:144
b0
static constexpr float b0
Definition: L1EGammaCrystalsEmulatorProducer.cc:80

◆ Pol2()

Double_t Pol2 ( Double_t *  x,
Double_t *  par 
)

Definition at line 1 of file fitMaximumPoint.cc.

References x.

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
float x
Definition: beamSpotDipStandalone.cc:55

◆ Pol2_Special()

Double_t Pol2_Special ( Double_t *  x,
Double_t *  par 
)

Definition at line 11 of file fitMaximumPoint.cc.

References x.

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
float x
Definition: beamSpotDipStandalone.cc:55
Generated for CMSSW Reference Manual by doxygen 1.8.14