BOSS 7.0.8
BESIII Offline Software System
Loading...
Searching...
No Matches
IniCalib.cpp
Go to the documentation of this file.
1#include "include/IniCalib.h"
2#include "include/fun.h"
3#include <cmath>
4
5#include "TF1.h"
6#include "TStyle.h"
7
8IniCalib::IniCalib(){
9 cout << "Calibration type: IniCalib" << endl;
10}
11
12IniCalib::~IniCalib(){
13}
14
15void IniCalib::init(TObjArray* hlist, MdcCosGeom* pGeom){
16 m_pGeom = pGeom;
17 char hname[200];
18 m_fdcom = new TFolder("mCommon", "mCommon");
19 hlist->Add(m_fdcom);
20
21 m_fdTmap = new TFolder("mThitmap", "mThitmap");
22 hlist->Add(m_fdTmap);
23
24 m_fdTraw = new TFolder("mTraw", "mTraw");
25 hlist->Add(m_fdTraw);
26
27 m_fdTrawCel = new TFolder("mTrawCell", "mTrawCell");
28 hlist->Add(m_fdTrawCel);
29
30 m_fdQmap = new TFolder("mQhitmap", "mQhitmap");
31 hlist->Add(m_fdQmap);
32
33 m_fdQraw = new TFolder("mQraw", "mQraw");
34 hlist->Add(m_fdQraw);
35
36 m_fdQrawCel = new TFolder("mQrawCell", "mQrawCell");
37 hlist->Add(m_fdQrawCel);
38
39 m_hLayerHitmapT = new TH1F("mT_Hitmap_Layer", "", 43, -0.5, 42.5);
40 m_fdcom->Add(m_hLayerHitmapT);
41
42 m_hWireHitMapT = new TH1F("mT_Hitmap_Wire", "", 6796, -0.5, 6795.5);
43 m_fdcom->Add(m_hWireHitMapT);
44
45 m_hLayerHitmapQ = new TH1F("mQ_Hitmap_Layer", "", 43, -0.5, 42.5);
46 m_fdcom->Add(m_hLayerHitmapQ);
47
48 m_hWireHitMapQ = new TH1F("mQ_Hitmap_Wire", "", 6796, -0.5, 6795.5);
49 m_fdcom->Add(m_hWireHitMapQ);
50
51 m_hTesAll = new TH1F("mTesAll", "", 750, 0, 1500);
52 m_fdcom->Add(m_hTesAll);
53
54 m_hTesCal = new TH1F("mTesCal", "", 750, 0, 1500);
55 m_fdcom->Add(m_hTesCal);
56
57 m_hTesFlag = new TH1F("mTes_Flag", "", 300, -0.5, 299.5);
58 m_fdcom->Add(m_hTesFlag);
59
60 for(int lay=0; lay<NLAYER; lay++){
61 int ncel = pGeom->getLayer(lay)->getNcell();
62
63 sprintf(hname, "mT_hitmap_Lay%02d", lay);
64 m_hlaymapT[lay] = new TH1F(hname, "", ncel, -0.5, (float)ncel-0.5);
65 m_fdTmap -> Add(m_hlaymapT[lay]);
66
67 sprintf(hname, "mTDC_Lay%02d", lay);
68 m_htdc[lay] = new TH1F(hname, "", 800, 0, 20000);
69 m_fdTraw -> Add(m_htdc[lay]);
70
71 sprintf(hname, "mTraw_Lay%02d", lay);
72 m_htraw[lay] = new TH1F(hname, "", 500, 0, 1000);
73 m_fdTraw -> Add(m_htraw[lay]);
74
75 sprintf(hname, "mQ_hitmap_Lay%02d", lay);
76 m_hlaymapQ[lay] = new TH1F(hname, "", ncel, -0.5, (float)ncel-0.5);
77 m_fdQmap -> Add(m_hlaymapQ[lay]);
78
79 sprintf(hname, "mQraw_Lay%02d", lay);
80 m_hqraw[lay] = new TH1F(hname, "", 2000, 0, 4000);
81 m_fdQraw -> Add(m_hqraw[lay]);
82 }
83
84 for(int wir=0; wir<NWIRE; wir++){
85 int lay = m_pGeom -> getWire(wir) -> getLayerId();
86 int cel = m_pGeom -> getWire(wir) -> getCellId();
87
88 sprintf(hname, "mTraw_%02d_%03d_%04d", lay, cel, wir);
89 m_htrawCel[wir] = new TH1F(hname, "", 300, 0, 600);
90 m_fdTrawCel -> Add(m_htrawCel[wir]);
91
92 sprintf(hname, "mQraw_%02d_%03d_%04d", lay, cel, wir);
93 m_hqrawCel[wir] = new TH1F(hname, "", 2000, 0, 4000);
94 m_fdQrawCel -> Add(m_hqrawCel[wir]);
95 }
96}
97
98void IniCalib::mergeHist(TFile* fhist){
99 TFolder* fdcom = (TFolder*)fhist->Get("Common");
100 TFolder* fdTmap = (TFolder*)fhist->Get("Thitmap");
101 TFolder* fdTraw = (TFolder*)fhist->Get("Traw");
102 TFolder* fdTrawCel = (TFolder*)fhist->Get("TrawCell");
103 TFolder* fdQmap = (TFolder*)fhist->Get("Qhitmap");
104 TFolder* fdQraw = (TFolder*)fhist->Get("Qraw");
105 TFolder* fdQrawCel = (TFolder*)fhist->Get("QrawCell");
106
107 char hname[200];
108 TH1F* hist;
109 hist = (TH1F*)fdcom->FindObjectAny("T_Hitmap_Layer");
110 m_hLayerHitmapT->Add(hist);
111
112 hist = (TH1F*)fdcom->FindObjectAny("T_Hitmap_Wire");
113 m_hWireHitMapT->Add(hist);
114
115 hist = (TH1F*)fdcom->FindObjectAny("Q_Hitmap_Layer");
116 m_hLayerHitmapQ->Add(hist);
117
118 hist = (TH1F*)fdcom->FindObjectAny("Q_Hitmap_Wire");
119 m_hWireHitMapQ->Add(hist);
120
121 hist = (TH1F*)fdcom->FindObjectAny("TesAll");
122 m_hTesAll->Add(hist);
123
124 hist = (TH1F*)fdcom->FindObjectAny("TesCal");
125 m_hTesCal->Add(hist);
126
127 hist = (TH1F*)fdcom->FindObjectAny("Tes_Flag");
128 m_hTesFlag->Add(hist);
129
130 for(int lay=0; lay<NLAYER; lay++){
131 sprintf(hname, "T_hitmap_Lay%02d", lay);
132 hist = (TH1F*)fdTmap->FindObjectAny(hname);
133 m_hlaymapT[lay]->Add(hist);
134
135 sprintf(hname, "TDC_Lay%02d", lay);
136 hist = (TH1F*)fdTraw->FindObjectAny(hname);
137 m_htdc[lay]->Add(hist);
138
139 sprintf(hname, "Traw_Lay%02d", lay);
140 hist = (TH1F*)fdTraw->FindObjectAny(hname);
141 m_htraw[lay]->Add(hist);
142
143 sprintf(hname, "Q_hitmap_Lay%02d", lay);
144 hist = (TH1F*)fdQmap->FindObjectAny(hname);
145 m_hlaymapQ[lay]->Add(hist);
146
147 sprintf(hname, "Qraw_Lay%02d", lay);
148 hist = (TH1F*)fdQraw->FindObjectAny(hname);
149 m_hqraw[lay]->Add(hist);
150 }
151
152 for(int wir=0; wir<NWIRE; wir++){
153 int lay = m_pGeom -> getWire(wir) -> getLayerId();
154 int cel = m_pGeom -> getWire(wir) -> getCellId();
155
156 sprintf(hname, "Traw_%02d_%03d_%04d", lay, cel, wir);
157 hist = (TH1F*)fdTrawCel->FindObjectAny(hname);
158 m_htrawCel[wir]->Add(hist);
159
160 sprintf(hname, "Qraw_%02d_%03d_%04d", lay, cel, wir);
161 hist = (TH1F*)fdQrawCel->FindObjectAny(hname);
162 m_hqrawCel[wir]->Add(hist);
163 }
164}
165
166void IniCalib::calib(MdcCalibConst* calconst, TObjArray* newXtList, TObjArray* r2tList){
167 int lay;
168 int wir;
169 double t0Fit[NLAYER];
170 double t0Cal[NLAYER];
171 double tmax[NLAYER];
172 double initT0 = gInitT0;
173
174 int fitTminFg[NLAYER];
175 int fitTmaxFg[NLAYER];
176 double chisq;
177 double ndf;
178 double chindfTmin[NLAYER];
179 double chindfTmax[NLAYER];
180 char funname[200];
181
182 // fit Tmin
183 TF1* ftmin[NLAYER];
184 for(lay=0; lay<NLAYER; lay++){
185 fitTminFg[lay] = 0;
186 chindfTmin[lay] = -1;
187 sprintf(funname, "ftmin%02d", lay);
188 ftmin[lay] = new TF1(funname, funTmin, 0, 150, 6);
189
190 if(1 == gFgCalib[lay]){
191 Stat_t nEntryTot = 0;
192 for(int ibin=1; ibin<=25; ibin++){
193 Stat_t entry = m_htraw[lay]->GetBinContent(ibin);
194 nEntryTot += entry;
195 }
196 double c0Ini = (double)nEntryTot / 25.0;
197 double c1Ini = (m_htraw[lay]->GetMaximum()) - c0Ini;
198
199 ftmin[lay] -> SetParameter(0, c0Ini);
200 ftmin[lay] -> SetParameter(1, c1Ini);
201 ftmin[lay] -> SetParameter(2, 0);
202 ftmin[lay] -> SetParameter(4, initT0);
203 ftmin[lay] -> SetParameter(5, 1);
204
205 m_htraw[lay] -> Fit(funname, "Q", "", gTminFitRange[lay][0], gTminFitRange[lay][1]);
206 gStyle -> SetOptFit(11);
207 chisq = ftmin[lay]->GetChisquare();
208 ndf = ftmin[lay]->GetNDF();
209 chindfTmin[lay] = chisq / ndf;
210 if(chindfTmin[lay] < gTminFitChindf){
211 fitTminFg[lay] = 1;
212 t0Fit[lay] = ftmin[lay]->GetParameter(4);
213 t0Fit[lay] += gT0Shift;
214 t0Cal[lay] = t0Fit[lay] - gTimeShift;
215 }
216 }
217
218 if(0 == fitTminFg[lay]){
219 wir = m_pGeom->getWire(lay, 0)->getWireId();
220 t0Cal[lay] = calconst->getT0(wir);
221 t0Fit[lay] = t0Cal[lay] + gTimeShift;
222 }
223 }
224
225 // fit Tmax
226 TF1* ftmax[NLAYER];
227 for(lay=0; lay<NLAYER; lay++){
228 fitTmaxFg[lay] = 0;
229 chindfTmax[lay] = -1;
230 sprintf(funname, "ftmax%02d", lay);
231 ftmax[lay] = new TF1(funname, funTmax, 250, 500, 4);
232
233 if(1 == gFgCalib[lay]){
234 ftmax[lay] -> SetParameter(2, gInitTm[lay]);
235 ftmax[lay] -> SetParameter(3, 10);
236
237 m_htraw[lay] -> Fit(funname, "Q+", "", gTmaxFitRange[lay][0], gTmaxFitRange[lay][1]);
238 gStyle -> SetOptFit(11);
239 chisq = ftmax[lay]->GetChisquare();
240 ndf = ftmax[lay]->GetNDF();
241 chindfTmax[lay] = chisq / ndf;
242 if(chindfTmax[lay] < gTmaxFitChindf){
243 fitTmaxFg[lay] = 1;
244 tmax[lay] = ftmax[lay]->GetParameter(2);
245 }
246 }
247
248 if(0 == fitTmaxFg[lay]){
249 tmax[lay] = (calconst->getXtpar(lay, 0, 0, 6)) + t0Fit[lay];
250 }
251 }
252
253 // output for check
254 ofstream ft0("iniT0.dat");
255 for(lay=0; lay<NLAYER; lay++){
256 ft0 << setw(5) << lay << setw(3) << fitTminFg[lay]
257 << setw(12) << t0Cal[lay] << setw(12) << t0Fit[lay]
258 << setw(12) << chindfTmin[lay] << setw(5) << fitTmaxFg[lay]
259 << setw(12) << tmax[lay] << setw(12) << tmax[lay] - t0Fit[lay]
260 << setw(12) << chindfTmax[lay] << endl;
261 }
262 ft0.close();
263 cout << "iniT0.dat was written." << endl;
264
265 // set T0
266 int i;
267 int nwire = m_pGeom -> getWireSize();
268 for(i=0; i<nwire; i++){
269 lay = m_pGeom -> getWire(i) -> getLayerId();
270 if(1 == gFgCalib[lay]){
271 calconst -> resetT0(i, t0Cal[lay]);
272 calconst -> resetDelT0(i, 0.0);
273 }
274 }
275
276 if(0 == gFgIniCalConst){
277 // set X-T
278 int lr;
279 int iEntr;
280 int ord;
281 double xtpar;
282 double xtini[8] = {0, 0.03, 0, 0, 0, 0, 999.9, 0};
283 for(lay=0; lay<NLAYER; lay++){
284 if(1 != gFgCalib[lay]) continue;
285
286 for(iEntr=0; iEntr<NENTRXT; iEntr++){
287 for(lr=0; lr<NLR; lr++){
288 for(ord=0; ord<NXTPAR; ord++){
289 if(6 == ord){
290 xtpar = tmax[lay] - t0Fit[lay];
291 } else{
292 xtpar = xtini[ord];
293 }
294 calconst -> resetXtpar(lay, iEntr, lr, ord, xtpar);
295 }
296 }
297 }
298 }
299
300 // set Q-T
301 for(lay=0; lay<NLAYER; lay++){
302 calconst -> resetQtpar0(lay, 0.0);
303 calconst -> resetQtpar1(lay, 0.0);
304 }
305
306 // set S-D
307 int bin;
308 double sdpar = 0.18; // mm
309 for(lay=0; lay<NLAYER; lay++){
310 for(iEntr=0; iEntr<NENTRSD; iEntr++){
311 for(lr=0; lr<2; lr++){
312 for(bin=0; bin<NSDBIN; bin++){
313 calconst -> resetSdpar(lay, iEntr, lr, bin, sdpar);
314 }
315 }
316 }
317 }
318 } else if(2 == gFgIniCalConst){
319 int lr;
320 int iEntr;
321 double xtpar;
322 for(lay=0; lay<NLAYER; lay++){
323 for(iEntr=0; iEntr<NENTRXT; iEntr++){
324 for(lr=0; lr<NLR; lr++){
325 xtpar = tmax[lay] - t0Fit[lay];
326 calconst -> resetXtpar(lay, iEntr, lr, 6, xtpar);
327 }
328 }
329 }
330 }
331
332 renameHist();
333 for(lay=0; lay<NLAYER; lay++){
334 delete ftmin[lay];
335 delete ftmax[lay];
336 }
337}
338
339Double_t IniCalib::funTmin(Double_t* x, Double_t* par){
340 Double_t fitval;
341 fitval = par[0] + par[1]*exp( -par[2]*(x[0]-par[3]) ) /
342 ( 1 + exp( -(x[0]-par[4])/par[5] ));
343 return fitval;
344}
345
346Double_t IniCalib::funTmax(Double_t* x, Double_t* par){
347 Double_t fitval;
348 fitval = par[0] + par[1] / (1 + exp((x[0]-par[2])/par[3]));
349 return fitval;
350}
351
352void IniCalib::renameHist(){
353 char hname[200];
354 m_fdcom->SetName("Common");
355 m_fdTmap->SetName("Thitmap");
356 m_fdTraw->SetName("Traw");
357 m_fdTrawCel->SetName("TrawCell");
358 m_fdQmap->SetName("Qhitmap");
359 m_fdQraw->SetName("Qraw");
360 m_fdQrawCel->SetName("QrawCell");
361
362 m_hLayerHitmapT->SetName("T_Hitmap_Layer");
363 m_hWireHitMapT->SetName("T_Hitmap_Wire");
364 m_hLayerHitmapQ->SetName("Q_Hitmap_Layer");
365 m_hWireHitMapQ->SetName("Q_Hitmap_Wire");
366 m_hTesAll->SetName("TesAll");
367 m_hTesCal->SetName("TesCal");
368 m_hTesFlag->SetName("Tes_Flag");
369
370 for(int lay=0; lay<NLAYER; lay++){
371 sprintf(hname, "T_hitmap_Lay%02d", lay);
372 m_hlaymapT[lay]->SetName(hname);
373
374 sprintf(hname, "TDC_Lay%02d", lay);
375 m_htdc[lay]->SetName(hname);
376
377 sprintf(hname, "Traw_Lay%02d", lay);
378 m_htraw[lay]->SetName(hname);
379
380 sprintf(hname, "Q_hitmap_Lay%02d", lay);
381 m_hlaymapQ[lay]->SetName(hname);
382
383 sprintf(hname, "Qraw_Lay%02d", lay);
384 m_hqraw[lay]->SetName(hname);
385 }
386 for(int wir=0; wir<NWIRE; wir++){
387 int lay = m_pGeom -> getWire(wir) -> getLayerId();
388 int cel = m_pGeom -> getWire(wir) -> getCellId();
389
390 sprintf(hname, "Traw_%02d_%03d_%04d", lay, cel, wir);
391 m_htrawCel[wir]->SetName(hname);
392
393 sprintf(hname, "Qraw_%02d_%03d_%04d", lay, cel, wir);
394 m_hqrawCel[wir]->SetName(hname);
395 }
396}
sprintf
sprintf(cut,"kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
Add
mg Add(gr3)
Fit
void Fit()
Definition: Eangle1D/Fit.cxx:3
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
IniCalib::calib
void calib(MdcCalibConst *calconst, TObjArray *newXtList, TObjArray *r2tList)
Definition: IniCalib.cpp:166
IniCalib::mergeHist
void mergeHist(TFile *fhist)
Definition: IniCalib.cpp:98
IniCalib::IniCalib
IniCalib()
Definition: IniCalib.cpp:8
IniCalib::init
void init(TObjArray *hlist, MdcCosGeom *pGeom)
Definition: IniCalib.cpp:15
IniCalib::~IniCalib
~IniCalib()
Definition: IniCalib.cpp:12
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