00001 #include "include/PreT0Calib.h"
00002 #include "include/fun.h"
00003 #include <cmath>
00004 #include "TF1.h"
00005 #include "TStyle.h"
00006
00007 PreT0Calib::PreT0Calib(){
00008 cout << "Calibration type: PreT0Calib" << endl;
00009 m_nzbin = 11;
00010 }
00011
00012 PreT0Calib::~PreT0Calib(){
00013 }
00014
00015 void PreT0Calib::init(TObjArray* hlist, MdcCosGeom* pGeom){
00016 m_pGeom = pGeom;
00017 char hname[200];
00018
00019 m_fdTrec = new TFolder("mTrec", "Trec");
00020 hlist->Add(m_fdTrec);
00021
00022 m_fdTrecZ = new TFolder("mTrecZ", "TrecZ");
00023 hlist->Add(m_fdTrecZ);
00024
00025 for(int lay=0; lay<NLAYER; lay++){
00026 for(int lr=0; lr<NLR; lr++){
00027 if(0 == lr) sprintf(hname, "mTrec%02d_L", lay);
00028 else if(1 == lr) sprintf(hname, "mTrec%02d_R", lay);
00029 else sprintf(hname, "mTrec%02d_C", lay);
00030
00031 if(lay < 8) m_hTrec[lay][lr] = new TH1F(hname, "", 100, 0, 400);
00032 else m_hTrec[lay][lr] = new TH1F(hname, "", 125, 0, 500);
00033 m_fdTrec -> Add(m_hTrec[lay][lr]);
00034 }
00035 }
00036
00037 for(int lay=0; lay<NLAYER; lay++){
00038 for(int iud=0; iud<2; iud++){
00039 if(0 == iud) sprintf(hname, "mTrecCosm%02d_Up", lay);
00040 else sprintf(hname, "mTrecCosm%02d_Dw", lay);
00041 if(lay < 8) m_hTrecCosm[lay][iud] = new TH1F(hname, "", 100, 0, 400);
00042 else m_hTrecCosm[lay][iud] = new TH1F(hname, "", 125, 0, 500);
00043 m_fdTrec -> Add(m_hTrecCosm[lay][iud]);
00044 }
00045 }
00046
00047 for(int lay=0; lay<NLAYER; lay++){
00048 for(int lr=0; lr<NLR; lr++){
00049 for(int bin=0; bin<m_nzbin; bin++){
00050 if(0 == lr) sprintf(hname, "mTrec%02d_z%02d_L", lay, bin);
00051 else if(1 == lr) sprintf(hname, "mTrec%02d_z%02d_R", lay, bin);
00052 else sprintf(hname, "mTrec%02d_z%02d_C", lay, bin);
00053
00054 if(lay < 8) m_hTrecZ[lay][lr][bin] = new TH1F(hname, "", 100, 0, 400);
00055 else m_hTrecZ[lay][lr][bin] = new TH1F(hname, "", 125, 0, 500);
00056 m_fdTrecZ -> Add(m_hTrecZ[lay][lr][bin]);
00057 }
00058 }
00059 }
00060
00061 }
00062
00063 void PreT0Calib::mergeHist(TFile* fhist){
00064 char hname[200];
00065 TFolder* fdTrec = (TFolder*)fhist->Get("Trec");
00066 TFolder* fdTrecZ = (TFolder*)fhist->Get("TrecZ");
00067
00068 TH1F* hist;
00069 for(int lay=0; lay<NLAYER; lay++){
00070 for(int lr=0; lr<NLR; lr++){
00071 if(0 == lr) sprintf(hname, "Trec%02d_L", lay);
00072 else if(1 == lr) sprintf(hname, "Trec%02d_R", lay);
00073 else sprintf(hname, "Trec%02d_C", lay);
00074 hist = (TH1F*)fdTrec->FindObjectAny(hname);
00075 m_hTrec[lay][lr]->Add(hist);
00076 }
00077 }
00078
00079 for(int lay=0; lay<NLAYER; lay++){
00080 for(int iud=0; iud<2; iud++){
00081 if(0 == iud) sprintf(hname, "TrecCosm%02d_Up", lay);
00082 else sprintf(hname, "TrecCosm%02d_Dw", lay);
00083 hist = (TH1F*)fdTrec->FindObjectAny(hname);
00084 m_hTrecCosm[lay][iud]->Add(hist);
00085 }
00086 }
00087
00088 for(int lay=0; lay<NLAYER; lay++){
00089 for(int lr=0; lr<NLR; lr++){
00090 for(int bin=0; bin<m_nzbin; bin++){
00091 if(0 == lr) sprintf(hname, "Trec%02d_z%02d_L", lay, bin);
00092 else if(1 == lr) sprintf(hname, "Trec%02d_z%02d_R", lay, bin);
00093 else sprintf(hname, "Trec%02d_z%02d_C", lay, bin);
00094 hist = (TH1F*)fdTrecZ->FindObjectAny(hname);
00095 m_hTrecZ[lay][lr][bin]->Add(hist);
00096 }
00097 }
00098 }
00099 }
00100
00101 void PreT0Calib::calib(MdcCalibConst* calconst, TObjArray* newXtList, TObjArray* r2tList){
00102
00103 double t0FitPar[NLAYER][NLR][6];
00104 double t0Fit[NLAYER][NLR];
00105 double t0Cal[NLAYER][NLR];
00106 double tmax[NLAYER][NLR];
00107 double initT0 = gInitT0;
00108
00109 int fitTminFg[NLAYER][NLR];
00110 int fitTmaxFg[NLAYER][NLR];
00111 double chisq;
00112 double ndf;
00113 double chindfTmin[NLAYER][NLR];
00114 double chindfTmax[NLAYER][NLR];
00115 char funname[200];
00116
00117
00118 TF1* ftminCosm[NLAYER][2];
00119 double t0FitCosm[NLAYER][2];
00120
00121 bool fgT0Ini = false;
00122 double t0ParIni[NLAYER][NLR][6];
00123 ifstream fparIni("fitT0_inival.txt");
00124 if(fparIni.is_open()){
00125 string strtmp;
00126 for(int lay=0; lay<NLAYER; lay++){
00127 fparIni >> strtmp >> strtmp;
00128 for(int i=0; i<6; i++) fparIni >> t0ParIni[lay][2][i];
00129 }
00130 fparIni.close();
00131 fgT0Ini = true;
00132 cout << "read initial values of T0 fit from fitT0_inival.txt" << endl;
00133 }
00134 if(!fgT0Ini) cout << "set initial values of T0 fit to default values" << endl;
00135
00136 TF1* ftmin[NLAYER][NLR];
00137 for(int lay=0; lay<NLAYER; lay++){
00138 for(int lr=0; lr<NLR; lr++){
00139 fitTminFg[lay][lr] = 0;
00140 chindfTmin[lay][lr] = -1;
00141 sprintf(funname, "ftmin%02d_%d", lay, lr);
00142 ftmin[lay][lr] = new TF1(funname, funTmin, 0, 150, 6);
00143 if(lr<2) continue;
00144
00145 if(1 == gFgCalib[lay]){
00146
00147
00148
00149
00150
00151
00152 double c1Ini = (m_hTrec[lay][lr]->GetMaximum());
00153
00154 if(fgT0Ini){
00155 for(int i=0; i<6; i++){
00156 if(fabs(t0ParIni[lay][2][i] + 9999)<0.01) continue;
00157 ftmin[lay][lr] -> SetParameter(i, t0ParIni[lay][2][i]);
00158 }
00159 } else{
00160 ftmin[lay][lr] -> SetParameter(0, 0);
00161 ftmin[lay][lr] -> SetParameter(1, c1Ini);
00162 ftmin[lay][lr] -> SetParameter(2, 0);
00163 ftmin[lay][lr] -> SetParameter(4, initT0);
00164 if(lay<4) ftmin[lay][lr] -> SetParameter(5, 4);
00165 else ftmin[lay][lr] -> SetParameter(5, 1.5);
00166 }
00167
00168 if(lay<4) m_hTrec[lay][lr] -> Fit(funname, "Q", "", 0, 140);
00169 else m_hTrec[lay][lr] -> Fit(funname, "Q", "", gTminFitRange[lay][0], gTminFitRange[lay][1]);
00170 gStyle -> SetOptFit(11);
00171 chisq = ftmin[lay][lr]->GetChisquare();
00172 ndf = ftmin[lay][lr]->GetNDF();
00173 chindfTmin[lay][lr] = chisq / ndf;
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 if(chindfTmin[lay][lr] < gTminFitChindf){
00184 fitTminFg[lay][lr] = 1;
00185 t0Fit[lay][lr] = ftmin[lay][lr]->GetParameter(4);
00186 t0Fit[lay][lr] += gT0Shift;
00187 t0Cal[lay][lr] = t0Fit[lay][lr] - gTimeShift;
00188 for(int i=0; i<6; i++) t0FitPar[lay][lr][i] = ftmin[lay][lr]->GetParameter(i);
00189 }
00190 }
00191
00192 if(0 == fitTminFg[lay][lr]){
00193 int wir = m_pGeom->getWire(lay, 0)->getWireId();
00194 t0Cal[lay][lr] = calconst->getT0(wir);
00195 t0Fit[lay][lr] = t0Cal[lay][lr] + gTimeShift;
00196 }
00197 }
00198
00199 for(int iud=0; iud<2; iud++){
00200 sprintf(funname, "ftminCosm_%02d_%d", lay, iud);
00201 ftminCosm[lay][iud] = new TF1(funname, funTmin, 0, 150, 6);
00202 ftminCosm[lay][iud] -> SetParameter(0, 0);
00203 ftminCosm[lay][iud] -> SetParameter(4, initT0);
00204 ftminCosm[lay][iud] -> SetParameter(5, 1);
00205 m_hTrecCosm[lay][iud] -> Fit(funname, "Q", "", gTminFitRange[lay][0], gTminFitRange[lay][1]);
00206 gStyle -> SetOptFit(11);
00207 t0FitCosm[lay][iud] += gT0Shift;
00208 t0FitCosm[lay][iud] = ftminCosm[lay][iud]->GetParameter(4);
00209 }
00210 }
00211
00212
00213 TF1* ftmax[NLAYER][NLR];
00214 for(int lay=0; lay<NLAYER; lay++){
00215 for(int lr=0; lr<NLR; lr++){
00216 fitTmaxFg[lay][lr] = 0;
00217 chindfTmax[lay][lr] = -1;
00218 sprintf(funname, "ftmax%02d_%d", lay, lr);
00219 ftmax[lay][lr] = new TF1(funname, funTmax, 250, 500, 4);
00220
00221 if(1 == gFgCalib[lay]){
00222 ftmax[lay][lr] -> SetParameter(2, gInitTm[lay]);
00223 ftmax[lay][lr] -> SetParameter(3, 10);
00224 m_hTrec[lay][lr] -> Fit(funname, "Q+", "", gTmaxFitRange[lay][0], gTmaxFitRange[lay][1]);
00225 gStyle -> SetOptFit(11);
00226 chisq = ftmax[lay][lr]->GetChisquare();
00227 ndf = ftmax[lay][lr]->GetNDF();
00228 chindfTmax[lay][lr] = chisq / ndf;
00229 if(chindfTmax[lay][lr] < gTmaxFitChindf){
00230 fitTmaxFg[lay][lr] = 1;
00231 tmax[lay][lr] = ftmax[lay][lr]->GetParameter(2);
00232 }
00233 }
00234
00235 if(0 == fitTmaxFg[lay][lr]){
00236 tmax[lay][lr] = (calconst->getXtpar(lay, 0, lr, 6)) + t0Fit[lay][2];
00237 }
00238 }
00239 }
00240
00241
00242 ofstream ft0("preT0.dat");
00243 for(int lay=0; lay<NLAYER; lay++){
00244 ft0 << setw(5) << lay << setw(3) << fitTminFg[lay][2]
00245 << setw(15) << t0Cal[lay][2] << setw(15) << t0Fit[lay][2]
00246 << setw(15) << chindfTmin[lay][2] << endl;
00247 }
00248 ft0 << endl;
00249 for(int lay=0; lay<NLAYER; lay++){
00250 ft0 << setw(5) << lay
00251 << setw(3) << fitTmaxFg[lay][0] << setw(10) << tmax[lay][0]
00252 << setw(10) << chindfTmax[lay][0]
00253 << setw(3) << fitTmaxFg[lay][1] << setw(10) << tmax[lay][1]
00254 << setw(10) << chindfTmax[lay][1]
00255 << setw(3) << fitTmaxFg[lay][2] << setw(10) << tmax[lay][2]
00256 << setw(10) << chindfTmax[lay][2]
00257 << setw(10) << tmax[lay][0] - t0Fit[lay][2]
00258 << setw(10) << tmax[lay][1] - t0Fit[lay][2]
00259 << setw(10) << tmax[lay][2] - t0Fit[lay][2]
00260 << endl;
00261 }
00262 ft0 << endl;
00263 for(int lay=0; lay<NLAYER; lay++){
00264 ft0 << setw(5) << lay << setw(15) << chindfTmin[lay][2];
00265 for(int i=0; i<6; i++) ft0 << setw(15) << t0FitPar[lay][2][i];
00266 ft0 << endl;
00267 }
00268 ft0.close();
00269 cout << "preT0.dat was written." << endl;
00270
00271
00272 ofstream ft0cosm("cosmicT0.dat");
00273 for(int lay=0; lay<NLAYER; lay++){
00274 ft0cosm << setw(5) << lay << setw(15) << t0Fit[lay][2]
00275 << setw(15) << t0FitCosm[lay][0] << setw(15) << t0FitCosm[lay][1] << endl;
00276 }
00277 ft0cosm.close();
00278
00279
00280 for(int i=0; i<NWIRE; i++){
00281 int lay = m_pGeom -> getWire(i) -> getLayerId();
00282 if(1 == gFgCalib[lay]){
00283 calconst -> resetT0(i, t0Cal[lay][2]);
00284 calconst -> resetDelT0(i, 0.0);
00285 }
00286 }
00287
00288
00289 if(gPreT0SetTm){
00290 double tm;
00291 for(int lay=0; lay<NLAYER; lay++){
00292 if(1 != gFgCalib[lay]) continue;
00293
00294 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
00295 for(int lr=0; lr<NLR; lr++){
00296 tm = tmax[lay][lr] - t0Fit[lay][2];
00297 if( (tmax[lay][lr] > gTmaxFitRange[lay][0]) &&
00298 (tmax[lay][lr] < gTmaxFitRange[lay][1]) ){
00299 calconst -> resetXtpar(lay, iEntr, lr, 6, tm);
00300 }
00301 }
00302 }
00303 }
00304 }
00305
00306 renameHist();
00307 for(int lay=0; lay<NLAYER; lay++){
00308 for(int lr=0; lr<NLR; lr++){
00309 delete ftmin[lay][lr];
00310 delete ftmax[lay][lr];
00311 }
00312 }
00313 }
00314
00315 void PreT0Calib::renameHist(){
00316 char hname[200];
00317 for(int lay=0; lay<NLAYER; lay++){
00318 for(int lr=0; lr<NLR; lr++){
00319 if(0 == lr) sprintf(hname, "Trec%02d_L", lay);
00320 else if(1 == lr) sprintf(hname, "Trec%02d_R", lay);
00321 else sprintf(hname, "Trec%02d_C", lay);
00322 m_hTrec[lay][lr]->SetName(hname);
00323 }
00324 }
00325 for(int lay=0; lay<NLAYER; lay++){
00326 for(int iud=0; iud<2; iud++){
00327 if(0 == iud) sprintf(hname, "TrecCosm%02d_Up", lay);
00328 else sprintf(hname, "TrecCosm%02d_Dw", lay);
00329 m_hTrecCosm[lay][iud]->SetName(hname);
00330 }
00331 }
00332 for(int lay=0; lay<NLAYER; lay++){
00333 for(int lr=0; lr<NLR; lr++){
00334 for(int bin=0; bin<m_nzbin; bin++){
00335 if(0 == lr) sprintf(hname, "Trec%02d_z%02d_L", lay, bin);
00336 else if(1 == lr) sprintf(hname, "Trec%02d_z%02d_R", lay, bin);
00337 else sprintf(hname, "Trec%02d_z%02d_C", lay, bin);
00338 m_hTrecZ[lay][lr][bin]->SetName(hname);
00339 }
00340 }
00341 }
00342 }
00343
00344 Double_t PreT0Calib::funTmin(Double_t* x, Double_t* par){
00345 Double_t fitval;
00346 fitval = par[0] + par[1]*exp( -par[2]*(x[0]-par[3]) ) /
00347 ( 1 + exp( -(x[0]-par[4])/par[5] ));
00348 return fitval;
00349 }
00350
00351 Double_t PreT0Calib::funTmax(Double_t* x, Double_t* par){
00352 Double_t fitval;
00353 fitval = par[0] + par[1] / (1 + exp((x[0]-par[2])/par[3]));
00354 return fitval;
00355 }