00001 #include "include/IniCalib.h"
00002 #include "include/fun.h"
00003 #include <cmath>
00004
00005 #include "TF1.h"
00006 #include "TStyle.h"
00007
00008 IniCalib::IniCalib(){
00009 cout << "Calibration type: IniCalib" << endl;
00010 }
00011
00012 IniCalib::~IniCalib(){
00013 }
00014
00015 void IniCalib::init(TObjArray* hlist, MdcCosGeom* pGeom){
00016 m_pGeom = pGeom;
00017 char hname[200];
00018 m_fdcom = new TFolder("mCommon", "mCommon");
00019 hlist->Add(m_fdcom);
00020
00021 m_fdTmap = new TFolder("mThitmap", "mThitmap");
00022 hlist->Add(m_fdTmap);
00023
00024 m_fdTraw = new TFolder("mTraw", "mTraw");
00025 hlist->Add(m_fdTraw);
00026
00027 m_fdTrawCel = new TFolder("mTrawCell", "mTrawCell");
00028 hlist->Add(m_fdTrawCel);
00029
00030 m_fdQmap = new TFolder("mQhitmap", "mQhitmap");
00031 hlist->Add(m_fdQmap);
00032
00033 m_fdQraw = new TFolder("mQraw", "mQraw");
00034 hlist->Add(m_fdQraw);
00035
00036 m_fdQrawCel = new TFolder("mQrawCell", "mQrawCell");
00037 hlist->Add(m_fdQrawCel);
00038
00039 m_hLayerHitmapT = new TH1F("mT_Hitmap_Layer", "", 43, -0.5, 42.5);
00040 m_fdcom->Add(m_hLayerHitmapT);
00041
00042 m_hWireHitMapT = new TH1F("mT_Hitmap_Wire", "", 6796, -0.5, 6795.5);
00043 m_fdcom->Add(m_hWireHitMapT);
00044
00045 m_hLayerHitmapQ = new TH1F("mQ_Hitmap_Layer", "", 43, -0.5, 42.5);
00046 m_fdcom->Add(m_hLayerHitmapQ);
00047
00048 m_hWireHitMapQ = new TH1F("mQ_Hitmap_Wire", "", 6796, -0.5, 6795.5);
00049 m_fdcom->Add(m_hWireHitMapQ);
00050
00051 m_hTesAll = new TH1F("mTesAll", "", 750, 0, 1500);
00052 m_fdcom->Add(m_hTesAll);
00053
00054 m_hTesCal = new TH1F("mTesCal", "", 750, 0, 1500);
00055 m_fdcom->Add(m_hTesCal);
00056
00057 m_hTesFlag = new TH1F("mTes_Flag", "", 300, -0.5, 299.5);
00058 m_fdcom->Add(m_hTesFlag);
00059
00060 for(int lay=0; lay<NLAYER; lay++){
00061 int ncel = pGeom->getLayer(lay)->getNcell();
00062
00063 sprintf(hname, "mT_hitmap_Lay%02d", lay);
00064 m_hlaymapT[lay] = new TH1F(hname, "", ncel, -0.5, (float)ncel-0.5);
00065 m_fdTmap -> Add(m_hlaymapT[lay]);
00066
00067 sprintf(hname, "mTDC_Lay%02d", lay);
00068 m_htdc[lay] = new TH1F(hname, "", 800, 0, 20000);
00069 m_fdTraw -> Add(m_htdc[lay]);
00070
00071 sprintf(hname, "mTraw_Lay%02d", lay);
00072 m_htraw[lay] = new TH1F(hname, "", 500, 0, 1000);
00073 m_fdTraw -> Add(m_htraw[lay]);
00074
00075 sprintf(hname, "mQ_hitmap_Lay%02d", lay);
00076 m_hlaymapQ[lay] = new TH1F(hname, "", ncel, -0.5, (float)ncel-0.5);
00077 m_fdQmap -> Add(m_hlaymapQ[lay]);
00078
00079 sprintf(hname, "mQraw_Lay%02d", lay);
00080 m_hqraw[lay] = new TH1F(hname, "", 2000, 0, 4000);
00081 m_fdQraw -> Add(m_hqraw[lay]);
00082 }
00083
00084 for(int wir=0; wir<NWIRE; wir++){
00085 int lay = m_pGeom -> getWire(wir) -> getLayerId();
00086 int cel = m_pGeom -> getWire(wir) -> getCellId();
00087
00088 sprintf(hname, "mTraw_%02d_%03d_%04d", lay, cel, wir);
00089 m_htrawCel[wir] = new TH1F(hname, "", 300, 0, 600);
00090 m_fdTrawCel -> Add(m_htrawCel[wir]);
00091
00092 sprintf(hname, "mQraw_%02d_%03d_%04d", lay, cel, wir);
00093 m_hqrawCel[wir] = new TH1F(hname, "", 2000, 0, 4000);
00094 m_fdQrawCel -> Add(m_hqrawCel[wir]);
00095 }
00096 }
00097
00098 void IniCalib::mergeHist(TFile* fhist){
00099 TFolder* fdcom = (TFolder*)fhist->Get("Common");
00100 TFolder* fdTmap = (TFolder*)fhist->Get("Thitmap");
00101 TFolder* fdTraw = (TFolder*)fhist->Get("Traw");
00102 TFolder* fdTrawCel = (TFolder*)fhist->Get("TrawCell");
00103 TFolder* fdQmap = (TFolder*)fhist->Get("Qhitmap");
00104 TFolder* fdQraw = (TFolder*)fhist->Get("Qraw");
00105 TFolder* fdQrawCel = (TFolder*)fhist->Get("QrawCell");
00106
00107 char hname[200];
00108 TH1F* hist;
00109 hist = (TH1F*)fdcom->FindObjectAny("T_Hitmap_Layer");
00110 m_hLayerHitmapT->Add(hist);
00111
00112 hist = (TH1F*)fdcom->FindObjectAny("T_Hitmap_Wire");
00113 m_hWireHitMapT->Add(hist);
00114
00115 hist = (TH1F*)fdcom->FindObjectAny("Q_Hitmap_Layer");
00116 m_hLayerHitmapQ->Add(hist);
00117
00118 hist = (TH1F*)fdcom->FindObjectAny("Q_Hitmap_Wire");
00119 m_hWireHitMapQ->Add(hist);
00120
00121 hist = (TH1F*)fdcom->FindObjectAny("TesAll");
00122 m_hTesAll->Add(hist);
00123
00124 hist = (TH1F*)fdcom->FindObjectAny("TesCal");
00125 m_hTesCal->Add(hist);
00126
00127 hist = (TH1F*)fdcom->FindObjectAny("Tes_Flag");
00128 m_hTesFlag->Add(hist);
00129
00130 for(int lay=0; lay<NLAYER; lay++){
00131 sprintf(hname, "T_hitmap_Lay%02d", lay);
00132 hist = (TH1F*)fdTmap->FindObjectAny(hname);
00133 m_hlaymapT[lay]->Add(hist);
00134
00135 sprintf(hname, "TDC_Lay%02d", lay);
00136 hist = (TH1F*)fdTraw->FindObjectAny(hname);
00137 m_htdc[lay]->Add(hist);
00138
00139 sprintf(hname, "Traw_Lay%02d", lay);
00140 hist = (TH1F*)fdTraw->FindObjectAny(hname);
00141 m_htraw[lay]->Add(hist);
00142
00143 sprintf(hname, "Q_hitmap_Lay%02d", lay);
00144 hist = (TH1F*)fdQmap->FindObjectAny(hname);
00145 m_hlaymapQ[lay]->Add(hist);
00146
00147 sprintf(hname, "Qraw_Lay%02d", lay);
00148 hist = (TH1F*)fdQraw->FindObjectAny(hname);
00149 m_hqraw[lay]->Add(hist);
00150 }
00151
00152 for(int wir=0; wir<NWIRE; wir++){
00153 int lay = m_pGeom -> getWire(wir) -> getLayerId();
00154 int cel = m_pGeom -> getWire(wir) -> getCellId();
00155
00156 sprintf(hname, "Traw_%02d_%03d_%04d", lay, cel, wir);
00157 hist = (TH1F*)fdTrawCel->FindObjectAny(hname);
00158 m_htrawCel[wir]->Add(hist);
00159
00160 sprintf(hname, "Qraw_%02d_%03d_%04d", lay, cel, wir);
00161 hist = (TH1F*)fdQrawCel->FindObjectAny(hname);
00162 m_hqrawCel[wir]->Add(hist);
00163 }
00164 }
00165
00166 void IniCalib::calib(MdcCalibConst* calconst, TObjArray* newXtList, TObjArray* r2tList){
00167 int lay;
00168 int wir;
00169 double t0Fit[NLAYER];
00170 double t0Cal[NLAYER];
00171 double tmax[NLAYER];
00172 double initT0 = gInitT0;
00173
00174 int fitTminFg[NLAYER];
00175 int fitTmaxFg[NLAYER];
00176 double chisq;
00177 double ndf;
00178 double chindfTmin[NLAYER];
00179 double chindfTmax[NLAYER];
00180 char funname[200];
00181
00182
00183 TF1* ftmin[NLAYER];
00184 for(lay=0; lay<NLAYER; lay++){
00185 fitTminFg[lay] = 0;
00186 chindfTmin[lay] = -1;
00187 sprintf(funname, "ftmin%02d", lay);
00188 ftmin[lay] = new TF1(funname, funTmin, 0, 150, 6);
00189
00190 if(1 == gFgCalib[lay]){
00191 Stat_t nEntryTot = 0;
00192 for(int ibin=1; ibin<=25; ibin++){
00193 Stat_t entry = m_htraw[lay]->GetBinContent(ibin);
00194 nEntryTot += entry;
00195 }
00196 double c0Ini = (double)nEntryTot / 25.0;
00197 double c1Ini = (m_htraw[lay]->GetMaximum()) - c0Ini;
00198
00199 ftmin[lay] -> SetParameter(0, c0Ini);
00200 ftmin[lay] -> SetParameter(1, c1Ini);
00201 ftmin[lay] -> SetParameter(2, 0);
00202 ftmin[lay] -> SetParameter(4, initT0);
00203 ftmin[lay] -> SetParameter(5, 1);
00204
00205 m_htraw[lay] -> Fit(funname, "Q", "", gTminFitRange[lay][0], gTminFitRange[lay][1]);
00206 gStyle -> SetOptFit(11);
00207 chisq = ftmin[lay]->GetChisquare();
00208 ndf = ftmin[lay]->GetNDF();
00209 chindfTmin[lay] = chisq / ndf;
00210 if(chindfTmin[lay] < gTminFitChindf){
00211 fitTminFg[lay] = 1;
00212 t0Fit[lay] = ftmin[lay]->GetParameter(4);
00213 t0Fit[lay] += gT0Shift;
00214 t0Cal[lay] = t0Fit[lay] - gTimeShift;
00215 }
00216 }
00217
00218 if(0 == fitTminFg[lay]){
00219 wir = m_pGeom->getWire(lay, 0)->getWireId();
00220 t0Cal[lay] = calconst->getT0(wir);
00221 t0Fit[lay] = t0Cal[lay] + gTimeShift;
00222 }
00223 }
00224
00225
00226 TF1* ftmax[NLAYER];
00227 for(lay=0; lay<NLAYER; lay++){
00228 fitTmaxFg[lay] = 0;
00229 chindfTmax[lay] = -1;
00230 sprintf(funname, "ftmax%02d", lay);
00231 ftmax[lay] = new TF1(funname, funTmax, 250, 500, 4);
00232
00233 if(1 == gFgCalib[lay]){
00234 ftmax[lay] -> SetParameter(2, gInitTm[lay]);
00235 ftmax[lay] -> SetParameter(3, 10);
00236
00237 m_htraw[lay] -> Fit(funname, "Q+", "", gTmaxFitRange[lay][0], gTmaxFitRange[lay][1]);
00238 gStyle -> SetOptFit(11);
00239 chisq = ftmax[lay]->GetChisquare();
00240 ndf = ftmax[lay]->GetNDF();
00241 chindfTmax[lay] = chisq / ndf;
00242 if(chindfTmax[lay] < gTmaxFitChindf){
00243 fitTmaxFg[lay] = 1;
00244 tmax[lay] = ftmax[lay]->GetParameter(2);
00245 }
00246 }
00247
00248 if(0 == fitTmaxFg[lay]){
00249 tmax[lay] = (calconst->getXtpar(lay, 0, 0, 6)) + t0Fit[lay];
00250 }
00251 }
00252
00253
00254 ofstream ft0("iniT0.dat");
00255 for(lay=0; lay<NLAYER; lay++){
00256 ft0 << setw(5) << lay << setw(3) << fitTminFg[lay]
00257 << setw(12) << t0Cal[lay] << setw(12) << t0Fit[lay]
00258 << setw(12) << chindfTmin[lay] << setw(5) << fitTmaxFg[lay]
00259 << setw(12) << tmax[lay] << setw(12) << tmax[lay] - t0Fit[lay]
00260 << setw(12) << chindfTmax[lay] << endl;
00261 }
00262 ft0.close();
00263 cout << "iniT0.dat was written." << endl;
00264
00265
00266 int i;
00267 int nwire = m_pGeom -> getWireSize();
00268 for(i=0; i<nwire; i++){
00269 lay = m_pGeom -> getWire(i) -> getLayerId();
00270 if(1 == gFgCalib[lay]){
00271 calconst -> resetT0(i, t0Cal[lay]);
00272 calconst -> resetDelT0(i, 0.0);
00273 }
00274 }
00275
00276 if(0 == gFgIniCalConst){
00277
00278 int lr;
00279 int iEntr;
00280 int ord;
00281 double xtpar;
00282 double xtini[8] = {0, 0.03, 0, 0, 0, 0, 999.9, 0};
00283 for(lay=0; lay<NLAYER; lay++){
00284 if(1 != gFgCalib[lay]) continue;
00285
00286 for(iEntr=0; iEntr<NENTRXT; iEntr++){
00287 for(lr=0; lr<NLR; lr++){
00288 for(ord=0; ord<NXTPAR; ord++){
00289 if(6 == ord){
00290 xtpar = tmax[lay] - t0Fit[lay];
00291 } else{
00292 xtpar = xtini[ord];
00293 }
00294 calconst -> resetXtpar(lay, iEntr, lr, ord, xtpar);
00295 }
00296 }
00297 }
00298 }
00299
00300
00301 for(lay=0; lay<NLAYER; lay++){
00302 calconst -> resetQtpar0(lay, 0.0);
00303 calconst -> resetQtpar1(lay, 0.0);
00304 }
00305
00306
00307 int bin;
00308 double sdpar = 0.18;
00309 for(lay=0; lay<NLAYER; lay++){
00310 for(iEntr=0; iEntr<NENTRSD; iEntr++){
00311 for(lr=0; lr<2; lr++){
00312 for(bin=0; bin<NSDBIN; bin++){
00313 calconst -> resetSdpar(lay, iEntr, lr, bin, sdpar);
00314 }
00315 }
00316 }
00317 }
00318 } else if(2 == gFgIniCalConst){
00319 int lr;
00320 int iEntr;
00321 double xtpar;
00322 for(lay=0; lay<NLAYER; lay++){
00323 for(iEntr=0; iEntr<NENTRXT; iEntr++){
00324 for(lr=0; lr<NLR; lr++){
00325 xtpar = tmax[lay] - t0Fit[lay];
00326 calconst -> resetXtpar(lay, iEntr, lr, 6, xtpar);
00327 }
00328 }
00329 }
00330 }
00331
00332 renameHist();
00333 for(lay=0; lay<NLAYER; lay++){
00334 delete ftmin[lay];
00335 delete ftmax[lay];
00336 }
00337 }
00338
00339 Double_t IniCalib::funTmin(Double_t* x, Double_t* par){
00340 Double_t fitval;
00341 fitval = par[0] + par[1]*exp( -par[2]*(x[0]-par[3]) ) /
00342 ( 1 + exp( -(x[0]-par[4])/par[5] ));
00343 return fitval;
00344 }
00345
00346 Double_t IniCalib::funTmax(Double_t* x, Double_t* par){
00347 Double_t fitval;
00348 fitval = par[0] + par[1] / (1 + exp((x[0]-par[2])/par[3]));
00349 return fitval;
00350 }
00351
00352 void IniCalib::renameHist(){
00353 char hname[200];
00354 m_fdcom->SetName("Common");
00355 m_fdTmap->SetName("Thitmap");
00356 m_fdTraw->SetName("Traw");
00357 m_fdTrawCel->SetName("TrawCell");
00358 m_fdQmap->SetName("Qhitmap");
00359 m_fdQraw->SetName("Qraw");
00360 m_fdQrawCel->SetName("QrawCell");
00361
00362 m_hLayerHitmapT->SetName("T_Hitmap_Layer");
00363 m_hWireHitMapT->SetName("T_Hitmap_Wire");
00364 m_hLayerHitmapQ->SetName("Q_Hitmap_Layer");
00365 m_hWireHitMapQ->SetName("Q_Hitmap_Wire");
00366 m_hTesAll->SetName("TesAll");
00367 m_hTesCal->SetName("TesCal");
00368 m_hTesFlag->SetName("Tes_Flag");
00369
00370 for(int lay=0; lay<NLAYER; lay++){
00371 sprintf(hname, "T_hitmap_Lay%02d", lay);
00372 m_hlaymapT[lay]->SetName(hname);
00373
00374 sprintf(hname, "TDC_Lay%02d", lay);
00375 m_htdc[lay]->SetName(hname);
00376
00377 sprintf(hname, "Traw_Lay%02d", lay);
00378 m_htraw[lay]->SetName(hname);
00379
00380 sprintf(hname, "Q_hitmap_Lay%02d", lay);
00381 m_hlaymapQ[lay]->SetName(hname);
00382
00383 sprintf(hname, "Qraw_Lay%02d", lay);
00384 m_hqraw[lay]->SetName(hname);
00385 }
00386 for(int wir=0; wir<NWIRE; wir++){
00387 int lay = m_pGeom -> getWire(wir) -> getLayerId();
00388 int cel = m_pGeom -> getWire(wir) -> getCellId();
00389
00390 sprintf(hname, "Traw_%02d_%03d_%04d", lay, cel, wir);
00391 m_htrawCel[wir]->SetName(hname);
00392
00393 sprintf(hname, "Qraw_%02d_%03d_%04d", lay, cel, wir);
00394 m_hqrawCel[wir]->SetName(hname);
00395 }
00396 }
