00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #ifndef PMT_FULL_MODEL_FIT_H
00015 #define PMT_FULL_MODEL_FIT_H
00016
00017 #include <stdlib.h>
00018 #include <iostream>
00019 #include <TROOT.h>
00020 #include <TApplication.h>
00021 #include <TVirtualX.h>
00022 #include <TGFileDialog.h>
00023 #include <TPostScript.h>
00024 #include <TStyle.h>
00025 #include <TSystem.h>
00026 #include <TRootHelpDialog.h>
00027 #include <TMinuit.h>
00028 #include <TVirtualFitter.h>
00029 #include "TPaveStats.h"
00030 #include "TText.h"
00031
00032 #include <fstream>
00033 #include <iomanip>
00034 #include <string.h>
00035 #include <math.h>
00036 #include <TMath.h>
00037 #include "TF1.h"
00038 #include "TCanvas.h"
00039 #include <TH1F.h>
00040
00041 #include "Event/ReadoutPmtCrate.h"
00042 #include "StatisticsSvc/IStatisticsSvc.h"
00043 #include "DataSvc/ICableSvc.h"
00044 #include "DataSvc/ICalibDataSvc.h"
00045 #include "Context/ServiceMode.h"
00046
00047
00048 #include "pmt_full_model.h"
00049
00050 using namespace std;
00051
00052
00053
00054 #define FULL_PMT_RESPONSE_MODEL
00055
00056
00057 #define MAX_PE_NUMBER 20
00058
00060
00061 Double_t x[4096];
00062 Double_t y[4096];
00063 Float_t num[4096];
00064 Float_t amplitude[4096];
00065
00066
00067
00069
00070
00071
00072
00073 Double_t mu;
00074 Double_t Q0;
00075 Double_t Delta0;
00076 Double_t Q1;
00077 Double_t Delta1;
00078 Double_t w;
00079 Double_t A;
00080
00081
00083
00084 extern Int_t NON_ZERO_MAX_NUM;
00085 extern Double_t EVENT_TOTLE_NUM;
00086 extern Double_t N;
00087 extern Double_t xlow;
00088 extern Double_t xhigh;
00089 extern long NNN;
00090
00091
00093
00094
00095
00096
00097
00098 Double_t fitf(Double_t *x, Double_t *par)
00099 {
00100 Double_t RESULT=0;
00101 int Nmax=int(N);
00102
00103 NNN++;
00104
00105 if(x[0]<xlow||x[0]>xhigh)
00106 return 0;
00107 else
00108 {
00109
00110 for(Int_t i=0;i<=Nmax;i++)
00111 {
00112 #ifdef FULL_PMT_RESPONSE_MODEL
00113
00114 RESULT=RESULT+SER_n(i,x,par);
00115 #endif
00116 }
00117
00118 RESULT=RESULT*EVENT_TOTLE_NUM;
00119 }
00120
00121 if(NNN%50000==0)
00122 cout<<"in fitf_IV..............."<<RESULT<<" "<<NNN<<endl;
00123 return RESULT;
00124 }
00125
00128
00129
00130
00131
00135 int pmt_full_model_fit(TH1F *histo1, double *par, double *par_err)
00136 {
00137
00138 NNN=0;
00139
00140 EVENT_TOTLE_NUM=par_err[0];
00141 if(EVENT_TOTLE_NUM<=0)
00142 {
00143 cout<<"total event number less than 0"<<endl;
00144 return -1;
00145 }
00146
00147 for(int i=0;i<4096;i++)
00148 {
00149 x[i]=0;
00150 y[i]=0;
00151 num[i]=0;
00152 amplitude[i]=0;
00153 }
00154 Int_t max_nonzero_bin_X=0;
00155 Int_t min_nonzero_bin_X=0;
00156
00157
00158
00159
00160
00161
00162 for(Int_t i=0;i<4096;i++)
00163 {
00164 amplitude[i]=histo1->GetBinContent(i);
00165 num[i]=int(histo1->GetBinCenter(i)+1);
00166 if(amplitude[i]!=0)
00167 max_nonzero_bin_X=int(histo1->GetBinCenter(i)+1);
00168 x[i]=int(histo1->GetBinCenter(i)+1);
00169 y[i]=amplitude[i];
00170 }
00171
00172
00173
00174 for(int i=0;i<histo1->GetMaximumBin()-1;i++)
00175 {
00176 if(amplitude[i]==0)
00177 {
00178 min_nonzero_bin_X=int(histo1->GetBinCenter(i)+1);
00179 }
00180 else
00181 {
00182 break;
00183 }
00184 }
00185
00186
00187
00188
00189
00190
00191 if(par[0]<=0.01||par[0]>=5)
00192 {
00193 cout<<" average_mu="<<par[0]<<" not suitable"<<endl;
00194 return -1;
00195 }
00196 else if(par[0]<=0.5)
00197 {
00198 N=6;
00199 }
00200 else if(par[0]<=1)
00201 {
00202 N=12;
00203 }
00204 else
00205 {
00206 N=MAX_PE_NUMBER;
00207 }
00208
00209
00210 double high_range=500;
00211
00212 xlow=min_nonzero_bin_X-1;
00213 xhigh=double(max_nonzero_bin_X-1);
00214
00215 mu=par[0];
00216 double max_content_bin_X=histo1->GetBinCenter(histo1->GetMaximumBin());
00217 cout<<"the max_content_bin_X="<<max_content_bin_X<<endl;
00218
00219 double max_content_X=histo1->GetMaximum();
00220 double bin_step=50;
00221 for(int t=0;t<10&&(max_content_bin_X+high_range<=max_nonzero_bin_X);t++)
00222 {
00223 if(histo1->GetBinContent(int(max_content_bin_X+high_range))<=0.01*max_content_X)
00224 {
00225 break;
00226 }
00227 else
00228 {
00229 high_range=high_range+bin_step;
00230 }
00231 }
00232
00233 if(TMath::Abs(max_content_bin_X-histo1->GetMean())>100)
00234 {
00235 cout<<"max_content_bin_X="<<max_content_bin_X<<" larger than the mean="<<histo1->GetMean()<<" more than 100"<<endl;
00236 return -1;
00237 }
00238
00239
00240 if(xhigh>max_content_bin_X-1+high_range)
00241 xhigh=max_content_bin_X-1+high_range;
00242
00243 if(xlow>=max_content_bin_X-1)
00244 {
00245 cout<<"--------------+++++++++------------------------"<<endl;
00246 cout<<"xlow="<<xlow<<" larger than maxBin="<<max_content_bin_X-1<<endl;
00247 cout<<"--------------+++++++++------------------------"<<endl;
00248 return -1;
00249 }
00250
00251 {
00252 Q0=par[1];
00253 if(Q0<0||Q0>xhigh||TMath::IsNaN(Q0))
00254 Q0=0;
00255 Delta0=par[2];
00256 if(Delta0<0||TMath::IsNaN(Delta0))
00257 Delta0=Q0*0.02;
00258 Q1=par[3];
00259 if(Q1<=0||Q1>xhigh||TMath::IsNaN(Q1)||Q1>(max_content_bin_X-1)*2)
00260 {
00261 if(max_content_bin_X-1-Q0>0)
00262 {
00263 Q1=max_content_bin_X-1-Q0;
00264 }
00265 else
00266 {
00267 cout<<"--------------+++++++++------------------------"<<endl;
00268 cout<<"max_content_bin_X="<<max_content_bin_X<<" lower than Q0="<<Q0<<endl;
00269 cout<<"--------------+++++++++------------------------"<<endl;
00270 return -1;
00271 }
00272 }
00273 Delta1=par[4];
00274 if(Delta1<=0||TMath::IsNaN(Delta1))
00275 Delta1=Q1*0.3;
00276 w=par[5];
00277 if(w<=0||TMath::IsNaN(w))
00278 w=0.02;
00279 A=par[6];
00280 if(A<=0||TMath::IsNaN(A))
00281 A=2.5*1/(Q1);
00282 }
00283
00284 if(xlow<Q0+3*Delta0)
00285 {
00286 cout<<"xlow="<<xlow<<" larger than 3sigma("<<Delta0<<") of Q0("<<Q0<<")"<<endl;
00287 mu=mu-mu*(xlow-Q0)/Delta0*0.05;
00288 if(mu<0.5*par[0])
00289 {
00290 mu=0.5*par[0];
00291 }
00292 cout<<"the corrected mu="<<mu<<"(the raw mu="<<par[0]<<")"<<endl;
00293 }
00294
00295 if(xlow>=max_content_bin_X-1)
00296 {
00297 cout<<"--------------+++++++++------------------------"<<endl;
00298 cout<<"xlow="<<xlow<<" larger than maxBin="<<max_content_bin_X-1<<endl;
00299 cout<<"--------------+++++++++------------------------"<<endl;
00300 return -1;
00301 }
00302 if(xhigh>max_content_bin_X-1+high_range)
00303 {
00304 xhigh=max_content_bin_X-1+high_range;
00305 }
00306
00307 cout<<"the data mu ="<<par[0]<<endl;
00308 cout<<"the p.e. maximum for fitting is :"<<N<<endl;
00309 cout<<"the fitting range is:("<<xlow<<","<<xhigh<<")"<<endl;
00310 cout<<"mu="<<mu<<endl;
00311 cout<<"Q0="<<Q0<<endl;
00312 cout<<"Delta0="<<Delta0<<endl;
00313 cout<<"Q1="<<Q1<<endl;
00314 cout<<"Delta1="<<Delta1<<endl;
00315 cout<<"w="<<w<<endl;
00316 cout<<"A="<<A<<endl;
00317
00318
00319
00320
00321 cout<<"=================================fitting=========================="<<endl;
00322
00323 cout<<"=================================TF1=========================="<<endl;
00324 TF1 *myfunc=new TF1("myfunc",fitf,0,4096,7);
00325 histo1->SetAxisRange(xlow,xhigh);
00326 myfunc->SetRange(xlow,xhigh);
00327 myfunc->SetParameters(mu,Q0,Delta0,Q1,Delta1,w,A);
00328
00329
00330
00331 myfunc->SetParLimits(1,Q0+0.01,Q0+0.01);
00332 myfunc->SetParLimits(2,Delta0,Delta0);
00333
00334 histo1->Fit("myfunc","BNR");
00335
00336
00337
00338 myfunc->GetParameters(par);
00339 for(int i=0;i<7;i++)
00340 {
00341 par_err[i]=myfunc->GetParError(i);
00342 }
00343
00344
00345 cout<<" -----------------------Completed!-----------------------"<<endl;
00346
00347 return 0;
00348
00349 }
00350
00351
00352 #endif //over of PMT_FULL_MODEL_FIT_H
00353
00354