CGEM BOSS 6.6.5.h
BESIII Offline Software System
Loading...
Searching...
No Matches
PreT0Calib.cpp
Go to the documentation of this file.
1#include "include/PreT0Calib.h"
2#include "include/fun.h"
3#include <cmath>
4#include "TF1.h"
5#include "TStyle.h"
6
7PreT0Calib::PreT0Calib(){
8 cout << "Calibration type: PreT0Calib" << endl;
9 m_nzbin = 11;
10}
11
14
15void PreT0Calib::init(TObjArray* hlist, MdcCosGeom* pGeom){
16 m_pGeom = pGeom;
17 char hname[200];
18
19 m_fdTrec = new TFolder("mTrec", "Trec");
20 hlist->Add(m_fdTrec);
21
22 m_fdTrecZ = new TFolder("mTrecZ", "TrecZ");
23 hlist->Add(m_fdTrecZ);
24
25 for(int lay=0; lay<NLAYER; lay++){
26 for(int lr=0; lr<NLR; lr++){
27 if(0 == lr) sprintf(hname, "mTrec%02d_L", lay);
28 else if(1 == lr) sprintf(hname, "mTrec%02d_R", lay);
29 else sprintf(hname, "mTrec%02d_C", lay);
30
31 if(lay < 8) m_hTrec[lay][lr] = new TH1F(hname, "", 100, 0, 400);
32 else m_hTrec[lay][lr] = new TH1F(hname, "", 125, 0, 500);
33 m_fdTrec -> Add(m_hTrec[lay][lr]);
34 }
35 }
36
37 for(int lay=0; lay<NLAYER; lay++){
38 for(int iud=0; iud<2; iud++){
39 if(0 == iud) sprintf(hname, "mTrecCosm%02d_Up", lay);
40 else sprintf(hname, "mTrecCosm%02d_Dw", lay);
41 if(lay < 8) m_hTrecCosm[lay][iud] = new TH1F(hname, "", 100, 0, 400);
42 else m_hTrecCosm[lay][iud] = new TH1F(hname, "", 125, 0, 500);
43 m_fdTrec -> Add(m_hTrecCosm[lay][iud]);
44 }
45 }
46
47 for(int lay=0; lay<NLAYER; lay++){
48 for(int lr=0; lr<NLR; lr++){
49 for(int bin=0; bin<m_nzbin; bin++){
50 if(0 == lr) sprintf(hname, "mTrec%02d_z%02d_L", lay, bin);
51 else if(1 == lr) sprintf(hname, "mTrec%02d_z%02d_R", lay, bin);
52 else sprintf(hname, "mTrec%02d_z%02d_C", lay, bin);
53
54 if(lay < 8) m_hTrecZ[lay][lr][bin] = new TH1F(hname, "", 100, 0, 400);
55 else m_hTrecZ[lay][lr][bin] = new TH1F(hname, "", 125, 0, 500);
56 m_fdTrecZ -> Add(m_hTrecZ[lay][lr][bin]);
57 }
58 }
59 }
60
61}
62
63void PreT0Calib::mergeHist(TFile* fhist){
64 char hname[200];
65 TFolder* fdTrec = (TFolder*)fhist->Get("Trec");
66 TFolder* fdTrecZ = (TFolder*)fhist->Get("TrecZ");
67
68 TH1F* hist;
69 for(int lay=0; lay<NLAYER; lay++){
70 for(int lr=0; lr<NLR; lr++){
71 if(0 == lr) sprintf(hname, "Trec%02d_L", lay);
72 else if(1 == lr) sprintf(hname, "Trec%02d_R", lay);
73 else sprintf(hname, "Trec%02d_C", lay);
74 hist = (TH1F*)fdTrec->FindObjectAny(hname);
75 m_hTrec[lay][lr]->Add(hist);
76 }
77 }
78
79 for(int lay=0; lay<NLAYER; lay++){
80 for(int iud=0; iud<2; iud++){
81 if(0 == iud) sprintf(hname, "TrecCosm%02d_Up", lay);
82 else sprintf(hname, "TrecCosm%02d_Dw", lay);
83 hist = (TH1F*)fdTrec->FindObjectAny(hname);
84 m_hTrecCosm[lay][iud]->Add(hist);
85 }
86 }
87
88 for(int lay=0; lay<NLAYER; lay++){
89 for(int lr=0; lr<NLR; lr++){
90 for(int bin=0; bin<m_nzbin; bin++){
91 if(0 == lr) sprintf(hname, "Trec%02d_z%02d_L", lay, bin);
92 else if(1 == lr) sprintf(hname, "Trec%02d_z%02d_R", lay, bin);
93 else sprintf(hname, "Trec%02d_z%02d_C", lay, bin);
94 hist = (TH1F*)fdTrecZ->FindObjectAny(hname);
95 m_hTrecZ[lay][lr][bin]->Add(hist);
96 }
97 }
98 }
99}
100
101void PreT0Calib::calib(MdcCalibConst* calconst, TObjArray* newXtList, TObjArray* r2tList){
102 // fit Tmin int lay;
103 double t0FitPar[NLAYER][NLR][6];
104 double t0Fit[NLAYER][NLR];
105 double t0Cal[NLAYER][NLR];
106 double tmax[NLAYER][NLR];
107 double initT0 = gInitT0;
108
109 int fitTminFg[NLAYER][NLR];
110 int fitTmaxFg[NLAYER][NLR];
111 double chisq;
112 double ndf;
113 double chindfTmin[NLAYER][NLR];
114 double chindfTmax[NLAYER][NLR];
115 char funname[200];
116
117 // add for cosmic-ray
118 TF1* ftminCosm[NLAYER][2];
119 double t0FitCosm[NLAYER][2];
120
121 bool fgT0Ini = false;
122 double t0ParIni[NLAYER][NLR][6];
123 ifstream fparIni("fitT0_inival.txt");
124 if(fparIni.is_open()){
125 string strtmp;
126 for(int lay=0; lay<NLAYER; lay++){
127 fparIni >> strtmp >> strtmp;
128 for(int i=0; i<6; i++) fparIni >> t0ParIni[lay][2][i];
129 }
130 fparIni.close();
131 fgT0Ini = true;
132 cout << "read initial values of T0 fit from fitT0_inival.txt" << endl;
133 }
134 if(!fgT0Ini) cout << "set initial values of T0 fit to default values" << endl;
135
136 TF1* ftmin[NLAYER][NLR];
137 for(int lay=0; lay<NLAYER; lay++){
138 for(int lr=0; lr<NLR; lr++){
139 fitTminFg[lay][lr] = 0;
140 chindfTmin[lay][lr] = -1;
141 sprintf(funname, "ftmin%02d_%d", lay, lr);
142 ftmin[lay][lr] = new TF1(funname, funTmin, 0, 150, 6);
143 if(lr<2) continue;
144
145 if(1 == gFgCalib[lay]){
146// Stat_t nEntryTot = 0;
147// for(int ibin=1; ibin<=25; ibin++){
148// Stat_t entry = m_hTrec[lay][lr]->GetBinContent(ibin);
149// nEntryTot += entry;
150// }
151// double c0Ini = (double)nEntryTot / 25.0;
152 double c1Ini = (m_hTrec[lay][lr]->GetMaximum());
153
154 if(fgT0Ini){
155 for(int i=0; i<6; i++){
156 if(fabs(t0ParIni[lay][2][i] + 9999)<0.01) continue;
157 ftmin[lay][lr] -> SetParameter(i, t0ParIni[lay][2][i]);
158 }
159 } else{
160 ftmin[lay][lr] -> SetParameter(0, 0);
161 ftmin[lay][lr] -> SetParameter(1, c1Ini);
162 ftmin[lay][lr] -> SetParameter(2, 0);
163 ftmin[lay][lr] -> SetParameter(4, initT0);
164 if(lay<4) ftmin[lay][lr] -> SetParameter(5, 4);
165 else ftmin[lay][lr] -> SetParameter(5, 1.5);
166 }
167
168 if(lay<4) m_hTrec[lay][lr] -> Fit(funname, "Q", "", 0, 140);
169 else m_hTrec[lay][lr] -> Fit(funname, "Q", "", gTminFitRange[lay][0], gTminFitRange[lay][1]);
170 gStyle -> SetOptFit(11);
171 chisq = ftmin[lay][lr]->GetChisquare();
172 ndf = ftmin[lay][lr]->GetNDF();
173 chindfTmin[lay][lr] = chisq / ndf;
174
175// if(chindfTmin[lay][lr] > 100){
176// ftmin[lay][lr] -> SetParameter(1, c1Ini+2000);
177// m_hTrec[lay][lr] -> Fit(funname, "Q", "", gTminFitRange[lay][0], gTminFitRange[lay][1]);
178// chisq = ftmin[lay][lr]->GetChisquare();
179// ndf = ftmin[lay][lr]->GetNDF();
180// chindfTmin[lay][lr] = chisq / ndf;
181// }
182
183 if(chindfTmin[lay][lr] < gTminFitChindf){
184 fitTminFg[lay][lr] = 1;
185 t0Fit[lay][lr] = ftmin[lay][lr]->GetParameter(4);
186 t0Fit[lay][lr] += gT0Shift;
187 t0Cal[lay][lr] = t0Fit[lay][lr] - gTimeShift;
188 for(int i=0; i<6; i++) t0FitPar[lay][lr][i] = ftmin[lay][lr]->GetParameter(i);
189 }
190 }
191
192 if(0 == fitTminFg[lay][lr]){
193 int wir = m_pGeom->getWire(lay, 0)->getWireId();
194 t0Cal[lay][lr] = calconst->getT0(wir);
195 t0Fit[lay][lr] = t0Cal[lay][lr] + gTimeShift;
196 }
197 }
198
199 for(int iud=0; iud<2; iud++){
200 sprintf(funname, "ftminCosm_%02d_%d", lay, iud);
201 ftminCosm[lay][iud] = new TF1(funname, funTmin, 0, 150, 6);
202 ftminCosm[lay][iud] -> SetParameter(0, 0);
203 ftminCosm[lay][iud] -> SetParameter(4, initT0);
204 ftminCosm[lay][iud] -> SetParameter(5, 1);
205 m_hTrecCosm[lay][iud] -> Fit(funname, "Q", "", gTminFitRange[lay][0], gTminFitRange[lay][1]);
206 gStyle -> SetOptFit(11);
207 t0FitCosm[lay][iud] += gT0Shift;
208 t0FitCosm[lay][iud] = ftminCosm[lay][iud]->GetParameter(4);
209 }
210 }
211
212 // fit Tmax
213 TF1* ftmax[NLAYER][NLR];
214 for(int lay=0; lay<NLAYER; lay++){
215 for(int lr=0; lr<NLR; lr++){
216 fitTmaxFg[lay][lr] = 0;
217 chindfTmax[lay][lr] = -1;
218 sprintf(funname, "ftmax%02d_%d", lay, lr);
219 ftmax[lay][lr] = new TF1(funname, funTmax, 250, 500, 4);
220
221 if(1 == gFgCalib[lay]){
222 ftmax[lay][lr] -> SetParameter(2, gInitTm[lay]);
223 ftmax[lay][lr] -> SetParameter(3, 10);
224 m_hTrec[lay][lr] -> Fit(funname, "Q+", "", gTmaxFitRange[lay][0], gTmaxFitRange[lay][1]);
225 gStyle -> SetOptFit(11);
226 chisq = ftmax[lay][lr]->GetChisquare();
227 ndf = ftmax[lay][lr]->GetNDF();
228 chindfTmax[lay][lr] = chisq / ndf;
229 if(chindfTmax[lay][lr] < gTmaxFitChindf){
230 fitTmaxFg[lay][lr] = 1;
231 tmax[lay][lr] = ftmax[lay][lr]->GetParameter(2);
232 }
233 }
234
235 if(0 == fitTmaxFg[lay][lr]){
236 tmax[lay][lr] = (calconst->getXtpar(lay, 0, lr, 6)) + t0Fit[lay][2];
237 }
238 }
239 }
240
241 // output for check
242 ofstream ft0("preT0.dat");
243 for(int lay=0; lay<NLAYER; lay++){
244 ft0 << setw(5) << lay << setw(3) << fitTminFg[lay][2]
245 << setw(15) << t0Cal[lay][2] << setw(15) << t0Fit[lay][2]
246 << setw(15) << chindfTmin[lay][2] << endl;
247 }
248 ft0 << endl;
249 for(int lay=0; lay<NLAYER; lay++){
250 ft0 << setw(5) << lay
251 << setw(3) << fitTmaxFg[lay][0] << setw(10) << tmax[lay][0]
252 << setw(10) << chindfTmax[lay][0]
253 << setw(3) << fitTmaxFg[lay][1] << setw(10) << tmax[lay][1]
254 << setw(10) << chindfTmax[lay][1]
255 << setw(3) << fitTmaxFg[lay][2] << setw(10) << tmax[lay][2]
256 << setw(10) << chindfTmax[lay][2]
257 << setw(10) << tmax[lay][0] - t0Fit[lay][2]
258 << setw(10) << tmax[lay][1] - t0Fit[lay][2]
259 << setw(10) << tmax[lay][2] - t0Fit[lay][2]
260 << endl;
261 }
262 ft0 << endl;
263 for(int lay=0; lay<NLAYER; lay++){
264 ft0 << setw(5) << lay << setw(15) << chindfTmin[lay][2];
265 for(int i=0; i<6; i++) ft0 << setw(15) << t0FitPar[lay][2][i];
266 ft0 << endl;
267 }
268 ft0.close();
269 cout << "preT0.dat was written." << endl;
270
271 // output for cosmic T0
272 ofstream ft0cosm("cosmicT0.dat");
273 for(int lay=0; lay<NLAYER; lay++){
274 ft0cosm << setw(5) << lay << setw(15) << t0Fit[lay][2]
275 << setw(15) << t0FitCosm[lay][0] << setw(15) << t0FitCosm[lay][1] << endl;
276 }
277 ft0cosm.close();
278
279 // set T0
280 for(int i=0; i<NWIRE; i++){
281 int lay = m_pGeom -> getWire(i) -> getLayerId();
282 if(1 == gFgCalib[lay]){
283 calconst -> resetT0(i, t0Cal[lay][2]);
284 calconst -> resetDelT0(i, 0.0);
285 }
286 }
287
288 // set tm of X-T
289 if(gPreT0SetTm){
290 double tm;
291 for(int lay=0; lay<NLAYER; lay++){
292 if(1 != gFgCalib[lay]) continue;
293
294 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
295 for(int lr=0; lr<NLR; lr++){
296 tm = tmax[lay][lr] - t0Fit[lay][2];
297 if( (tmax[lay][lr] > gTmaxFitRange[lay][0]) &&
298 (tmax[lay][lr] < gTmaxFitRange[lay][1]) ){
299 calconst -> resetXtpar(lay, iEntr, lr, 6, tm);
300 }
301 }
302 }
303 }
304 }
305
306 renameHist();
307 for(int lay=0; lay<NLAYER; lay++){
308 for(int lr=0; lr<NLR; lr++){
309 delete ftmin[lay][lr];
310 delete ftmax[lay][lr];
311 }
312 }
313}
314
315void PreT0Calib::renameHist(){
316 char hname[200];
317 for(int lay=0; lay<NLAYER; lay++){
318 for(int lr=0; lr<NLR; lr++){
319 if(0 == lr) sprintf(hname, "Trec%02d_L", lay);
320 else if(1 == lr) sprintf(hname, "Trec%02d_R", lay);
321 else sprintf(hname, "Trec%02d_C", lay);
322 m_hTrec[lay][lr]->SetName(hname);
323 }
324 }
325 for(int lay=0; lay<NLAYER; lay++){
326 for(int iud=0; iud<2; iud++){
327 if(0 == iud) sprintf(hname, "TrecCosm%02d_Up", lay);
328 else sprintf(hname, "TrecCosm%02d_Dw", lay);
329 m_hTrecCosm[lay][iud]->SetName(hname);
330 }
331 }
332 for(int lay=0; lay<NLAYER; lay++){
333 for(int lr=0; lr<NLR; lr++){
334 for(int bin=0; bin<m_nzbin; bin++){
335 if(0 == lr) sprintf(hname, "Trec%02d_z%02d_L", lay, bin);
336 else if(1 == lr) sprintf(hname, "Trec%02d_z%02d_R", lay, bin);
337 else sprintf(hname, "Trec%02d_z%02d_C", lay, bin);
338 m_hTrecZ[lay][lr][bin]->SetName(hname);
339 }
340 }
341 }
342}
343
344Double_t PreT0Calib::funTmin(Double_t* x, Double_t* par){
345 Double_t fitval;
346 fitval = par[0] + par[1]*exp( -par[2]*(x[0]-par[3]) ) /
347 ( 1 + exp( -(x[0]-par[4])/par[5] ));
348 return fitval;
349}
350
351Double_t PreT0Calib::funTmax(Double_t* x, Double_t* par){
352 Double_t fitval;
353 fitval = par[0] + par[1] / (1 + exp((x[0]-par[2])/par[3]));
354 return fitval;
355}
Fit
gr Fit("g1","R")
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:252
bin
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition FoamA.h:85
MdcCalibConst::getT0
double getT0(int wireid) const
Definition MdcCalibAlg/MdcCalibConst.h:27
MdcCalibConst::getXtpar
double getXtpar(int lay, int entr, int lr, int order)
Definition MdcCalibConst.cpp:31
PreT0Calib::calib
void calib(MdcCalibConst *calconst, TObjArray *newXtList, TObjArray *r2tList)
Definition PreT0Calib.cpp:101
PreT0Calib::mergeHist
void mergeHist(TFile *fhist)
Definition PreT0Calib.cpp:63
PreT0Calib::init
void init(TObjArray *hlist, MdcCosGeom *pGeom)
Definition PreT0Calib.cpp:15