00001 #ifndef SaturationFcn_h 00002 #define SaturationFcn_h 00003 00004 #include "Minuit2/FCNBase.h" 00005 #include <vector> 00006 00007 class SaturationFcn : public FCNBase{ 00008 00009 public: 00010 00011 SaturationFcn(){} 00012 00013 ~SaturationFcn(){} 00014 00015 00016 void set_data(int N,float *charge_ptr,float *adc_ptr){ 00017 00018 float x[20],y[20]; 00019 00020 for(int i=0;i<N;i++){ 00021 x[i]=charge_ptr[i]; 00022 y[i]=adc_ptr[i]; 00023 datx[i]=x[i]; 00024 daty[i]=y[i]; 00025 //printf("%d datx daty %f %f \n",i,datx[i],daty[i]); 00026 } 00027 00028 x3start=(y[4]*x[1]-y[1]*x[4])/(x[1]-x[4]); 00029 x0start=daty[13]-x3start; 00030 x1start=(y[4]-y[1])/(x[4]-x[1])/x0start; 00031 x2start=20.; 00032 //printf(" x0-2start %f %f %f %f\n",x0start,x1start,x2start,x3start); 00033 } 00034 00035 virtual double Up() const {return 1.;} 00036 00037 virtual double operator()(const std::vector<double>& x) const { 00038 double chisq = 0.0; 00039 int N=20; 00040 for(int i=0;i<N;i++){ 00041 double val=1.0+pow(x[1]*datx[i],x[2]); 00042 double val2=1.0/x[2]; 00043 val=x[0]*x[1]*datx[i]/pow(val,val2); 00044 double tmp=(daty[i]-x[3]-val); 00045 //printf(" dat: %d %f %f %f %f \n",i,datx[i],daty[i]-x[3],val,tmp); 00046 chisq=chisq+tmp*tmp; 00047 } 00048 //printf("x0-3 %f %f %f %f chisq %f \n",x[0],x[1],x[2],x[3],chisq); 00049 return chisq; 00050 } 00051 00052 double x0start; 00053 double x1start; 00054 double x2start; 00055 double x3start; 00056 00057 private: 00058 double datx[20],daty[20]; 00059 }; 00060 00061 #endif