GFlash/TB/fitMaximumPoint.cc

0001 Double_t Pol2(Double_t *x, Double_t *par)
0002 {
0003   //polynoom(2)
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   //polynoom(2)
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){  // xxx
0025 //====================================================
0026 //====================================================
0027
0028   // -------------------------------------------
0029   // 1) Clean beam profile
0030   //     o all bins:  only bins >  2  entries
0031   //     o per x-bin: only bins > 75% of maximum
0032   // -------------------------------------------
0033   // get some caracteristics
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   // [1a] Only keep bins with more than 2 entries
0046   //      Also remember for each x-bin the maximum bin in y that satisfies this result
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   // [1b] Only keep events with more than 85% of the maximum
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   // 2) Find region to fit
0069   //     o Make profile by compuing means in every bin (error = RMS/sqrt(N))
0070   //     o Find maximum
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   // [2a] Make TH1F containing the mean values
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       //Int_t Nevt_slice = 0;
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   // [2b] Find maximum
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   // [2c] Define the fit range  (0.975% and error should be less than 8\%)
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   // 3) Do the actual fit
0144   //    o make clone of histogram
0145   // -------------------------------------------------
0146   TH1F *h1  = (TH1F *) IHBeamProf_clean_3->Clone();
0147   // 3a] Do the actual fit
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   //h1->Write(); //write the profile
0172   // ----------------------------
0173   // 4) Compute uncertainty on x0
0174   // ----------------------------
0175   // [4a] compute dxo using the covariance matrix
0176   // covariance matrix
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   // jacobiaan
0186   Float_t j1  =  b/(2*a*a);
0187   Float_t j2  = -1./(2*a);
0188   // determinant covariance martix
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   // inverse matrix
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   // variance
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   // 5) Second type of fit
0204   // ---------------------
0205   // 5a] Do the actual fit
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   //printf("Cov0[1][1] = %f\n",CovMat0[1][1]);
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 }