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