File indexing completed on 2023-03-17 11:24:57
0001 Double_t Pol2(Double_t *x, Double_t *par)
0002 {
0003
0004 Double_t part0 = par[0];
0005 Double_t part1 = par[1]*x[0];
0006 Double_t part2 = par[2]*x[0]*x[0];
0007 Double_t fitval = part0 + part1 + part2;
0008 return fitval;
0009 }
0010
0011 Double_t Pol2_Special(Double_t *x, Double_t *par)
0012 {
0013
0014 Double_t part0 = par[0];
0015 Double_t part1 = -2.*par[2]*par[1]*x[0];
0016 Double_t part2 = par[2]*x[0]*x[0];
0017 Double_t fitval = part0 + part1 + part2;
0018 return fitval;
0019 }
0020
0021
0022
0023
0024 Float_t Fit_MaximumPoint( TH2F *h2, Int_t Ireturn = 0){
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034 Int_t x_nbins = h2->GetXaxis()->GetNbins();
0035 Float_t x_max = h2->GetXaxis()->GetXmax();
0036 Float_t x_min = h2->GetXaxis()->GetXmin();
0037 Float_t delta_x = h2->GetXaxis()->GetBinWidth(10);
0038 Int_t y_nbins = h2->GetYaxis()->GetNbins();
0039 Float_t y_max = h2->GetYaxis()->GetXmax();
0040 Float_t y_min = h2->GetYaxis()->GetXmin();
0041 Float_t delta_y = h2->GetYaxis()->GetBinWidth(10);
0042 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);
0043 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);
0044
0045
0046
0047 Int_t y_binmax_high[500] = {1};
0048 for (Int_t ix=1;ix<x_nbins+1;ix++){
0049 for (Int_t iy=1;iy<y_nbins+1;iy++){
0050 Int_t Nevts_clean = 0;
0051 if( h2->GetBinContent(ix,iy) > 1 ) {
0052 IHBeamProf_clean_1->SetBinContent(ix,iy, h2->GetBinContent(ix,iy));
0053 if( iy > y_binmax_high[ix] ){ y_binmax_high[ix] = iy; }
0054 }
0055 }
0056 }
0057
0058 Int_t y_binmax_low[500] = {1};
0059 for (Int_t ix=1;ix<x_nbins+1;ix++){
0060 y_binmax_low[ix] = (Int_t)( 0.85*(Float_t)( y_binmax_high[ix]));
0061 for (Int_t iy=y_binmax_low[ix];iy<y_binmax_high[ix]+1;iy++){
0062 IHBeamProf_clean_2->SetBinContent(ix,iy,IHBeamProf_clean_1->GetBinContent(ix,iy));
0063 }
0064 }
0065
0066
0067
0068
0069
0070
0071
0072 TH1F *IHBeamProf_clean_3 = new TH1F((h2->GetName()+(TString("_clprof"))).Data(), "Meam Energy versus X (clean) ", x_nbins, x_min, x_max);
0073
0074
0075 TH1F *h_temp = new TH1F("h_temp"," Energy distribution for single bin in X", y_nbins, y_min, y_max);
0076 for (Int_t ix=1;ix<x_nbins+1;ix++){
0077 double Nevt_slice = 0;
0078
0079 for (Int_t iy=1;iy<y_nbins+1;iy++){
0080 Int_t Nevt_tmp = (int)IHBeamProf_clean_2->GetBinContent(ix,iy);
0081 Nevt_slice += Nevt_tmp;
0082 h_temp->SetBinContent(iy,Nevt_tmp);
0083 }
0084 Float_t Y_mean = h_temp->GetMean();
0085 Float_t Y_rms = h_temp->GetRMS();
0086 Float_t Y_mean_error = (Nevt_slice>0) ? Y_rms/sqrt(Nevt_slice) : 9999.;
0087 IHBeamProf_clean_3->SetBinContent(ix,Y_mean);
0088 IHBeamProf_clean_3->SetBinError(ix,Y_mean_error);
0089 printf("%d %f %f %d %f\n",ix,Y_mean,Y_rms,Nevt_slice,Y_mean_error);
0090 }
0091
0092
0093 Float_t Y_tmp_max = -999.;
0094 Float_t x_max_guess = 0.;
0095 Int_t x_max_guess_bin = 0;
0096 for (Int_t ix=1;ix<x_nbins+1;ix++){
0097 Float_t X_tmp = IHBeamProf_clean_3->GetBinCenter(ix);
0098 Float_t Y_tmp = IHBeamProf_clean_3->GetBinContent(ix);
0099 Float_t Y_tmp_error = IHBeamProf_clean_3->GetBinError(ix);
0100 if( (fabs(X_tmp) < 10.) && (Y_tmp > Y_tmp_max) && (Y_tmp_error < 100.)){
0101 Y_tmp_max = Y_tmp;
0102 x_max_guess = X_tmp;
0103 x_max_guess_bin = ix;
0104 printf("Xtmp = %5.3f Ytmp = %5.3f Ytmp_max = %5.3f\n",X_tmp,Y_tmp,Y_tmp_max);
0105 }
0106 }
0107 printf("Fit: o Start for maximum = %5.3f\n",x_max_guess);
0108
0109
0110 Int_t fitbinmin = 9999;
0111 Int_t fitbinmax = -9999;
0112 Float_t x_tmp,y_tmp,e_tmp,er_tmp;
0113 for (Int_t ix = x_max_guess_bin; ix<x_nbins+1;ix++){
0114 fitbinmax = ix;
0115 x_tmp = IHBeamProf_clean_3->GetBinCenter(ix);
0116 y_tmp = IHBeamProf_clean_3->GetBinContent(ix);
0117 e_tmp = IHBeamProf_clean_3->GetBinError(ix);
0118 er_tmp = (y_tmp>0.) ? (e_tmp/y_tmp) : 1.00;
0119 printf("%3d %f %f %f -- %f %f\n",ix,x_tmp,y_tmp,e_tmp,y_tmp/Y_tmp_max,er_tmp);
0120 if( y_tmp < 0.975*Y_tmp_max && er_tmp < 0.008) break;
0121 }
0122 for (Int_t ix = x_max_guess_bin; ix>1 ;ix--){
0123 fitbinmin = ix;
0124 x_tmp = IHBeamProf_clean_3->GetBinCenter(ix);
0125 y_tmp = IHBeamProf_clean_3->GetBinContent(ix);
0126 e_tmp = IHBeamProf_clean_3->GetBinError(ix);
0127 er_tmp = (y_tmp>0.) ? (e_tmp/y_tmp) : 1.00;
0128 printf("%3d %f %f %f -- %f %f\n",ix,x_tmp,y_tmp,e_tmp,y_tmp/Y_tmp_max,er_tmp);
0129 if( y_tmp < 0.975* Y_tmp_max && er_tmp < 0.008) break;
0130 }
0131
0132 Double_t posmin_fit = x_min+fitbinmin*delta_x;
0133 Double_t posmax_fit = x_min+fitbinmax*delta_x;
0134 printf(" o Fit range = %5.3f -- %5.3f -- %5.3f\n",posmin_fit,x_max_guess,posmax_fit);
0135 if( fabs(posmax_fit - posmin_fit ) < 4. ){
0136 printf("Something is wrong with this range: returning dummy value\n");
0137 posmin_fit = x_max_guess - 6.0;
0138 posmax_fit = x_max_guess + 6.0;
0139 return -99.00;
0140 }
0141
0142
0143
0144
0145
0146 TH1F *h1 = (TH1F *) IHBeamProf_clean_3->Clone();
0147
0148 Double_t fitresults[3] = {0.};
0149 Double_t fitresults_errors[3] = {0.};
0150 TF1 *f1 = new TF1("f1",Pol2,-50.,50.,3);
0151 f1->SetParameters( 1.,1.,1.);
0152 h1->Fit(f1,"Q","",posmin_fit, posmax_fit);
0153 for(int i=0 ; i< 3 ; i++) {
0154 fitresults[i] = f1->GetParameter(i);
0155 fitresults_errors[i] = f1->GetParError(i);
0156 }
0157 Float_t chi2 = f1->GetChisquare()/f1->GetNDF();
0158 Float_t a = fitresults[2]; Float_t da = fitresults_errors[2];
0159 Float_t b = fitresults[1]; Float_t db = fitresults_errors[1];
0160 Float_t c = fitresults[0]; Float_t dc = fitresults_errors[0];
0161 Float_t x0 = (-1.*b)/(2*a);
0162 printf(" a = %7.2f b = %7.2f c = %7.2f\n",a,b,c);
0163 printf(" da = %7.2f db = %7.2f dc = %7.2f\n",da,db,dc);
0164
0165 cout << "risultati del fit polinomiale: " << fitresults[0] << " " << fitresults[1] << " " << fitresults[2] << endl;
0166
0167 char myProfTitle[200];
0168 sprintf(myProfTitle, h2->GetTitle());
0169 strcat(myProfTitle,"_prof");
0170 h1->Write(myProfTitle);
0171
0172
0173
0174
0175
0176
0177 Double_t CovMat[3][3];
0178 gMinuit->mnemat(&CovMat[0][0],3);
0179 Float_t v11 = CovMat[0][0]; Float_t v12 = CovMat[0][1]; Float_t v13 = CovMat[0][2];
0180 Float_t v21 = CovMat[1][0]; Float_t v22 = CovMat[1][1]; Float_t v23 = CovMat[1][2];
0181 Float_t v31 = CovMat[2][0]; Float_t v32 = CovMat[2][1]; Float_t v33 = CovMat[2][2];
0182 printf("Covariance Matrix: v11 = %f v12 = %f v13 = %f\n",CovMat[0][0],CovMat[0][1],CovMat[0][2]);
0183 printf(" v21 = %f v22 = %f v23 = %f\n",CovMat[1][0],CovMat[1][1],CovMat[1][2]);
0184 printf(" v31 = %f v32 = %f v33 = %f\n",CovMat[2][0],CovMat[2][1],CovMat[2][2]);
0185
0186 Float_t j1 = b/(2*a*a);
0187 Float_t j2 = -1./(2*a);
0188
0189 Float_t det = v11*(v33*v22-v32*v23)-v21*(v33*v12-v32*v13)+v31*(v23*v12-v22*v13);
0190 printf("Determinant = %f\n",det);
0191
0192 Float_t v11_i = v33*v22-v32*v23; Float_t v12_i = -(v33*v12-v32*v13); Float_t v13_i = v23*v12-v22*v13;
0193 Float_t v21_i = -(v33*v21-v31*v23); Float_t v22_i = v33*v11-v31*v13 ; Float_t v23_i = -(v23*v11-v21*v13) ;
0194 Float_t v31_i = v32*v21-v31*v22; Float_t v32_i = -(v32*v11-v31*v12); Float_t v33_i = v22*v11-v21*v12;
0195
0196 Float_t var = j1*(j1*v11_i+j2*v12_i)+j2*(j1*v21_i+j2*v22_i);
0197 var /= det;
0198 Float_t dx0 = sqrt(var);
0199 printf("Type 1 fit: o x0 = %f +/- %f\n",x0,dx0);
0200
0201
0202
0203
0204
0205
0206 TH1F *h0 = (TH1F *) IHBeamProf_clean_3->Clone();
0207 Double_t fitresults0[3] = {0.};
0208 Double_t fitresults0_errors[3] = {0.};
0209 TF1 *f0 = new TF1("f0",Pol2_Special,-50.,50.,3);
0210 f0->SetParameters( -1.,x_max_guess,5000.);
0211 h0->Fit(f0,"Q","",posmin_fit, posmax_fit);
0212 for(int i=0 ; i< 3 ; i++) {
0213 fitresults0[i] = f0->GetParameter(i);
0214 fitresults0_errors[i] = f0->GetParError(i);
0215 }
0216 Float_t a0 = fitresults0[2]; Float_t da0 = fitresults0_errors[2];
0217 Float_t b0 = fitresults0[1]; Float_t db0 = fitresults0_errors[1];
0218 Float_t c0 = fitresults0[0]; Float_t dc0 = fitresults0_errors[0];
0219
0220 Double_t CovMat0[3][3];
0221 gMinuit->mnemat(&CovMat0[0][0],3);
0222
0223 Float_t db0cov = sqrt(CovMat0[1][1]);
0224 printf("Type 2 fit: o x0 = %f +/- %f %f \n",b0,db0,db0cov);
0225 delete IHBeamProf_clean_1;
0226 delete IHBeamProf_clean_2;
0227 delete IHBeamProf_clean_3;
0228 delete h_temp;
0229 delete f1;
0230 delete f0;
0231
0232 if(Ireturn == 10) {return x0;}
0233 if(Ireturn == 11) {return dx0;}
0234 if(Ireturn == 20) {return b0;}
0235 if(Ireturn == 21) {return db0cov;}
0236 if(Ireturn == 30) {return chi2;}
0237
0238 if(Ireturn == 99) {return c0;}
0239
0240
0241
0242 }