BOSS 7.0.6
BESIII Offline Software System
Loading...
Searching...
No Matches
PreT0MdcCalib Class Reference

#include <PreT0MdcCalib.h>

+ Inheritance diagram for PreT0MdcCalib:

Public Member Functions

 PreT0MdcCalib ()
 
 ~PreT0MdcCalib ()
 
void initialize (TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)
 
void setParam (MdcCalParams &param)
 
int fillHist (MdcCalEvent *event)
 
int updateConst (MdcCalibConst *calconst)
 
void clear ()
 
- Public Member Functions inherited from MdcCalib
 MdcCalib ()
 
virtual ~MdcCalib ()
 
virtual void initialize (TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)=0
 
virtual void setParam (MdcCalParams &param)=0
 
virtual int fillHist (MdcCalEvent *event)=0
 
virtual int updateConst (MdcCalibConst *calconst)=0
 
virtual void clear ()=0
 

Detailed Description

Definition at line 14 of file PreT0MdcCalib.h.

Constructor & Destructor Documentation

◆ PreT0MdcCalib()

PreT0MdcCalib::PreT0MdcCalib ( )

Definition at line 31 of file PreT0MdcCalib.cxx.

31 {
32 m_nzbin = 11;
33}

◆ ~PreT0MdcCalib()

PreT0MdcCalib::~PreT0MdcCalib ( )

Definition at line 35 of file PreT0MdcCalib.cxx.

35 {
36}

Member Function Documentation

◆ clear()

void PreT0MdcCalib::clear ( )
virtual

Implements MdcCalib.

Definition at line 38 of file PreT0MdcCalib.cxx.

38 {
39 int lay;
40 int lr;
41 int bin;
42 for(lay=0; lay<MdcCalNLayer; lay++){
43 for(lr=0; lr<MdcCalLR; lr++){
44 delete m_hTrec[lay][lr];
45 for(bin=0; bin<m_nzbin; bin++){
46 delete m_hTrecZ[lay][lr][bin];
47 }
48 }
49 }
50 for(lay=0; lay<MdcCalNLayer; lay++){
51 for(bin=0; bin<2; bin++) delete m_hTrecCosm[lay][bin];
52 }
53 delete m_fdTrec;
54 delete m_fdTrecZ;
55
57}
*******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
const int MdcCalNLayer
Definition: MdcCalParams.h:6
const int MdcCalLR
Definition: MdcCalParams.h:10
virtual void clear()=0
Definition: MdcCalib.cxx:76

◆ fillHist()

int PreT0MdcCalib::fillHist ( MdcCalEvent event)
virtual

Implements MdcCalib.

Definition at line 138 of file PreT0MdcCalib.cxx.

138 {
139 IMessageSvc* msgSvc;
140 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
141 MsgStream log(msgSvc, "PreT0MdcCalib");
142 log << MSG::DEBUG << "PreT0MdcCalib::fillHist()" << endreq;
143
144// MdcCalib::fillHist(event);
145
146 int i;
147 int k;
148 int lay;
149 int cel;
150 int lr;
151 int bin;
152
153 double tdc;
154 double traw;
155 double trec;
156 double tdr;
157 double t0;
158 double tp = 0.0;
159 double tof = 0.0;
160 double zhit;
161 double zprop;
162
163 MdcCalRecTrk* rectrk;
164 MdcCalRecHit* rechit;
165
166 IDataProviderSvc* eventSvc = NULL;
167 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
168
169 // get EsTimeCol
170 double tes = event->getTes();
171 // cout << "tes " << tes << endl;
172 bool esCutFg = event->getEsCutFlag();
173 if( ! esCutFg ) return -1;
174
175 int nhit;
176 int ntrk = event -> getNTrk();
177 for(i=0; i<ntrk; i++){
178 rectrk = event->getRecTrk(i);
179 nhit = rectrk -> getNHits();
180 for(k=0; k<nhit; k++){
181 rechit = rectrk -> getRecHit(k);
182 lay = rechit -> getLayid();
183 cel = rechit -> getCellid();
184 lr = rechit -> getLR();
185 zhit = rechit -> getZhit();
186 tdc = rechit -> getTdc();
187 tdr = rechit -> getTdrift();
188
189 traw = tdc * MdcCalTdcCnv;
190// traw = tdc; // ns
191
192 zprop = fabs(zhit - m_zst[lay]);
193 bin = (int)(zprop / m_zwid[lay]);
194 tp = zprop / m_vp[lay];
195 t0 = m_mdcFunSvc -> getT0(lay, cel);
196 tof = rechit -> getTof();
197
198 double adc = rechit -> getQhit();
199 double tw = m_mdcFunSvc->getTimeWalk(lay, adc);
200// cout << setw(15) << adc << setw(15) << tw << endl;
201
202 const MdcGeoWire* pWire = m_mdcGeomSvc -> Wire(lay, cel);
203 double y1 = pWire -> Backward().y();
204 double y2 = pWire -> Forward().y();
205 double ymid = (y1+y2)*0.5;
206// if(5==m_param.particle){ // cosmic-ray
207// if(ymid > 0) trec = traw - tp + tof - tes + m_param.timeShift - tw;
208// else trec = traw - tp - tof - tes + m_param.timeShift - tw;
209// } else{
210// trec = traw - tp - tof - tes + m_param.timeShift - tw;
211// }
212 trec = traw - tp - tof - tes + m_param.timeShift - tw;
213
214 m_hTrec[lay][lr] -> Fill(trec);
215 m_hTrec[lay][2] -> Fill(trec); // l-r in common
216 if(ymid > 0) m_hTrecCosm[lay][0] -> Fill(trec);
217 else m_hTrecCosm[lay][1] -> Fill(trec);
218
219 if(bin < m_nzbin){
220 m_hTrecZ[lay][lr][bin] -> Fill(trec);
221 m_hTrecZ[lay][2][bin] -> Fill(trec);
222 }
223 }
224 }
225 return 1;
226}
curve Fill()
const double MdcCalTdcCnv
Definition: MdcCalParams.h:24
IMessageSvc * msgSvc()
virtual double getTimeWalk(int layid, double Q) const =0
double timeShift
Definition: MdcCalParams.h:42

◆ initialize()

void PreT0MdcCalib::initialize ( TObjArray *  hlist,
IMdcGeomSvc mdcGeomSvc,
IMdcCalibFunSvc mdcFunSvc,
IMdcUtilitySvc mdcUtilitySvc 
)
virtual

Implements MdcCalib.

Definition at line 59 of file PreT0MdcCalib.cxx.

60 {
61 IMessageSvc* msgSvc;
62 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
63 MsgStream log(msgSvc, "PreT0MdcCalib");
64 log << MSG::INFO << "PreT0MdcCalib::initialize()" << endreq;
65
66 m_hlist = hlist;
67 m_mdcGeomSvc = mdcGeomSvc;
68 m_mdcFunSvc = mdcFunSvc;
69 m_mdcUtilitySvc = mdcUtilitySvc;
70
71// MdcCalib::initialize(m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc);
72
73 int lay;
74 int lr;
75 int bin;
76 char hname[200];
77
78 m_fdTrec = new TFolder("Trec", "Trec");
79 m_hlist->Add(m_fdTrec);
80
81 m_fdTrecZ = new TFolder("TrecZ", "TrecZ");
82 m_hlist->Add(m_fdTrecZ);
83
84 for(lay=0; lay<MdcCalNLayer; lay++){
85 for(lr=0; lr<MdcCalLR; lr++){
86 if(0 == lr) sprintf(hname, "Trec%02d_L", lay);
87 else if(1 == lr) sprintf(hname, "Trec%02d_R", lay);
88 else sprintf(hname, "Trec%02d_C", lay);
89
90 if(lay < 8) m_hTrec[lay][lr] = new TH1F(hname, "", 100, 0, 400);
91 else m_hTrec[lay][lr] = new TH1F(hname, "", 125, 0, 500);
92 m_fdTrec -> Add(m_hTrec[lay][lr]);
93 }
94 }
95
96 for(lay=0; lay<MdcCalNLayer; lay++){
97 for(int iud=0; iud<2; iud++){
98 if(0 == iud) sprintf(hname, "TrecCosm%02d_Up", lay);
99 else sprintf(hname, "TrecCosm%02d_Dw", lay);
100 if(lay < 8) m_hTrecCosm[lay][iud] = new TH1F(hname, "", 100, 0, 400);
101 else m_hTrecCosm[lay][iud] = new TH1F(hname, "", 125, 0, 500);
102 m_fdTrec -> Add(m_hTrecCosm[lay][iud]);
103 }
104 }
105
106 for(lay=0; lay<MdcCalNLayer; lay++){
107 for(lr=0; lr<MdcCalLR; lr++){
108 for(bin=0; bin<m_nzbin; bin++){
109 if(0 == lr) sprintf(hname, "Trec%02d_z%02d_L", lay, bin);
110 else if(1 == lr) sprintf(hname, "Trec%02d_z%02d_R", lay, bin);
111 else sprintf(hname, "Trec%02d_z%02d_C", lay, bin);
112
113 if(lay < 8) m_hTrecZ[lay][lr][bin] = new TH1F(hname, "", 100, 0, 400);
114 else m_hTrecZ[lay][lr][bin] = new TH1F(hname, "", 125, 0, 500);
115 m_fdTrecZ -> Add(m_hTrecZ[lay][lr][bin]);
116 }
117 }
118 }
119
120 double zeast;
121 double zwest;
122 for(lay=0; lay<MdcCalNLayer; lay++){
123 zwest = m_mdcGeomSvc->Wire(lay, 0)->Forward().z();
124 zeast = m_mdcGeomSvc->Wire(lay, 0)->Backward().z();
125 m_zwid[lay] = (zeast - zwest) / (double)m_nzbin;
126
127 if(lay < 8) m_vp[lay] = 220.0; // *10^9 mm/s
128 else m_vp[lay] = 240.0; // *10^9 mm/s
129
130 if( 0 == (lay % 2) ){ // west end
131 m_zst[lay] = zwest;
132 } else{ // east end
133 m_zst[lay] = zeast;
134 }
135 }
136}
virtual const MdcGeoWire *const Wire(unsigned id)=0
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128
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)
mg Add(gr3)

◆ setParam()

void PreT0MdcCalib::setParam ( MdcCalParams param)
inlinevirtual

Implements MdcCalib.

Definition at line 51 of file PreT0MdcCalib.h.

51 {
52 MdcCalib::setParam(param);
53 m_param = param;
54}
virtual void setParam(MdcCalParams &param)=0
Definition: MdcCalib.h:293

◆ updateConst()

int PreT0MdcCalib::updateConst ( MdcCalibConst calconst)
virtual

Implements MdcCalib.

Definition at line 228 of file PreT0MdcCalib.cxx.

228 {
229 IMessageSvc* msgSvc;
230 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
231 MsgStream log(msgSvc, "PreT0MdcCalib");
232 log << MSG::DEBUG << "PreT0MdcCalib::updateConst()" << endreq;
233
234// MdcCalib::updateConst(calconst);
235
236 // fit Tmin
237 int lay;
238 int wir;
239 int lr;
240 double t0Fit[MdcCalNLayer][MdcCalLR];
241 double t0Cal[MdcCalNLayer][MdcCalLR];
242 double tmax[MdcCalNLayer][MdcCalLR];
243 double initT0 = m_param.initT0;
244
245 int fitTminFg[MdcCalNLayer][MdcCalLR];
246 int fitTmaxFg[MdcCalNLayer][MdcCalLR];
247 double chisq;
248 double ndf;
249 double chindfTmin[MdcCalNLayer][MdcCalLR];
250 double chindfTmax[MdcCalNLayer][MdcCalLR];
251 char funname[200];
252
253 // add for cosmic-ray
254 TF1* ftminCosm[MdcCalNLayer][2];
255 double t0FitCosm[MdcCalNLayer][2];
256
257 TF1* ftmin[MdcCalNLayer][MdcCalLR];
258// sprintf(funname, "ftmin");
259// TF1* ftmin = new TF1(funname, funTmin, 0, 150, 6);
260 for(lay=0; lay<MdcCalNLayer; lay++){
261 for(lr=0; lr<MdcCalLR; lr++){
262 fitTminFg[lay][lr] = 0;
263 chindfTmin[lay][lr] = -1;
264 sprintf(funname, "ftmin%02d_%d", lay, lr);
265 ftmin[lay][lr] = new TF1(funname, funTmin, 0, 150, 6);
266
267 if(1 == m_param.fgCalib[lay]){
268 Stat_t nEntryTot = 0;
269 for(int ibin=1; ibin<=25; ibin++){
270 Stat_t entry = m_hTrec[lay][lr]->GetBinContent(ibin);
271 nEntryTot += entry;
272 }
273 double c0Ini = (double)nEntryTot / 25.0;
274 double c1Ini = (m_hTrec[lay][lr]->GetMaximum()) - c0Ini;
275
276 ftmin[lay][lr] -> SetParameter(0, c0Ini);
277 ftmin[lay][lr] -> SetParameter(1, c1Ini);
278 ftmin[lay][lr] -> SetParameter(2, 0);
279 ftmin[lay][lr] -> SetParameter(4, initT0);
280 ftmin[lay][lr] -> SetParameter(5, 2);
281
282 m_hTrec[lay][lr] -> Fit(funname, "Q", "",
283 m_param.tminFitRange[lay][0],
284 m_param.tminFitRange[lay][1]);
285 gStyle -> SetOptFit(11);
286// chisq = ftmin[lay][lr]->GetChisquare();
287// ndf = ftmin[lay][lr]->GetNDF();
288 chisq = ftmin[lay][lr]->GetChisquare();
289 ndf = ftmin[lay][lr]->GetNDF();
290 chindfTmin[lay][lr] = chisq / ndf;
291
292 if(chindfTmin[lay][lr] < m_param.tminFitChindf){
293 fitTminFg[lay][lr] = 1;
294 t0Fit[lay][lr] = ftmin[lay][lr]->GetParameter(4);
295// t0Fit[lay][lr] = ftmin->GetParameter(4);
296 t0Fit[lay][lr] += m_param.t0Shift;
297 t0Cal[lay][lr] = t0Fit[lay][lr] - m_param.timeShift;
298 }
299 }
300
301 if(0 == fitTminFg[lay][lr]){
302 wir = m_mdcGeomSvc->Wire(lay, 0)->Id();
303 t0Cal[lay][lr] = calconst->getT0(wir);
304 t0Fit[lay][lr] = t0Cal[lay][lr] + m_param.timeShift;
305 }
306 }
307
308 for(int iud=0; iud<2; iud++){
309 sprintf(funname, "ftminCosm_%02d_%d", lay, iud);
310 ftminCosm[lay][iud] = new TF1(funname, funTmin, 0, 150, 6);
311 ftminCosm[lay][iud] -> SetParameter(0, 0);
312 ftminCosm[lay][iud] -> SetParameter(4, initT0);
313 ftminCosm[lay][iud] -> SetParameter(5, 1);
314 m_hTrecCosm[lay][iud] -> Fit(funname, "Q", "",
315 m_param.tminFitRange[lay][0],
316 m_param.tminFitRange[lay][1]);
317 gStyle -> SetOptFit(11);
318 t0FitCosm[lay][iud] += m_param.t0Shift;
319 t0FitCosm[lay][iud] = ftminCosm[lay][iud]->GetParameter(4);
320 }
321 }
322
323 // fit Tmax
324 TF1* ftmax[MdcCalNLayer][MdcCalLR];
325 for(lay=0; lay<MdcCalNLayer; lay++){
326 for(lr=0; lr<MdcCalLR; lr++){
327 fitTmaxFg[lay][lr] = 0;
328 chindfTmax[lay][lr] = -1;
329 sprintf(funname, "ftmax%02d_%d", lay, lr);
330 ftmax[lay][lr] = new TF1(funname, funTmax, 250, 500, 4);
331
332 if(1 == m_param.fgCalib[lay]){
333 ftmax[lay][lr] -> SetParameter(2, m_param.initTm[lay]);
334 ftmax[lay][lr] -> SetParameter(3, 10);
335 m_hTrec[lay][lr] -> Fit(funname, "Q+", "",
336 m_param.tmaxFitRange[lay][0],
337 m_param.tmaxFitRange[lay][1]);
338 gStyle -> SetOptFit(11);
339 chisq = ftmax[lay][lr]->GetChisquare();
340 ndf = ftmax[lay][lr]->GetNDF();
341 chindfTmax[lay][lr] = chisq / ndf;
342 if(chindfTmax[lay][lr] < m_param.tmaxFitChindf){
343 fitTmaxFg[lay][lr] = 1;
344 tmax[lay][lr] = ftmax[lay][lr]->GetParameter(2);
345 }
346 }
347
348 if(0 == fitTmaxFg[lay][lr]){
349 tmax[lay][lr] = (calconst->getXtpar(lay, 0, lr, 6)) + t0Fit[lay][2];
350 }
351 }
352 }
353
354 // output for check
355 ofstream ft0("preT0.dat");
356 for(lay=0; lay<MdcCalNLayer; lay++){
357 ft0 << setw(5) << lay << setw(3) << fitTminFg[lay][2]
358 << setw(15) << t0Cal[lay][2] << setw(15) << t0Fit[lay][2]
359 << setw(15) << chindfTmin[lay][2] << endl;
360 }
361 ft0 << endl;
362 for(lay=0; lay<MdcCalNLayer; lay++){
363 ft0 << setw(5) << lay
364 << setw(3) << fitTmaxFg[lay][0] << setw(10) << tmax[lay][0]
365 << setw(10) << chindfTmax[lay][0]
366 << setw(3) << fitTmaxFg[lay][1] << setw(10) << tmax[lay][1]
367 << setw(10) << chindfTmax[lay][1]
368 << setw(3) << fitTmaxFg[lay][2] << setw(10) << tmax[lay][2]
369 << setw(10) << chindfTmax[lay][2]
370 << setw(10) << tmax[lay][0] - t0Fit[lay][2]
371 << setw(10) << tmax[lay][1] - t0Fit[lay][2]
372 << setw(10) << tmax[lay][2] - t0Fit[lay][2]
373 << endl;
374 }
375 ft0.close();
376 cout << "preT0.dat was written." << endl;
377
378 // output for cosmic T0
379 ofstream ft0cosm("cosmicT0.dat");
380 for(lay=0; lay<MdcCalNLayer; lay++){
381 ft0cosm << setw(5) << lay << setw(15) << t0Fit[lay][2]
382 << setw(15) << t0FitCosm[lay][0] << setw(15) << t0FitCosm[lay][1] << endl;
383 }
384 ft0cosm.close();
385
386 // set T0
387 int i;
388 int nwire = m_mdcGeomSvc -> getWireSize();
389 for(i=0; i<nwire; i++){
390 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
391 if(1 == m_param.fgCalib[lay]){
392 calconst -> resetT0(i, t0Cal[lay][2]);
393 calconst -> resetDelT0(i, 0.0);
394 }
395 }
396
397 // set tm of X-T
398 if(m_param.preT0SetTm){
399 int iEntr;
400 double tm;
401 for(lay=0; lay<MdcCalNLayer; lay++){
402 if(1 != m_param.fgCalib[lay]) continue;
403
404 for(iEntr=0; iEntr<MdcCalNENTRXT; iEntr++){
405 for(lr=0; lr<MdcCalLR; lr++){
406 tm = tmax[lay][lr] - t0Fit[lay][2];
407 if( (tmax[lay][lr] > m_param.tmaxFitRange[lay][0]) &&
408 (tmax[lay][lr] < m_param.tmaxFitRange[lay][1]) ){
409 calconst -> resetXtpar(lay, iEntr, lr, 6, tm);
410 }
411 }
412 }
413 }
414 }
415
416 // set sigma
417 int bin;
418 double sdpar = m_param.initSigma; // mm
419 for(lay=0; lay<MdcCalNLayer; lay++){
420 for(int iEntr=0; iEntr<MdcCalNENTRSD; iEntr++){
421 for(lr=0; lr<2; lr++){
422 for(bin=0; bin<MdcCalSdNBIN; bin++){
423 calconst -> resetSdpar(lay, iEntr, lr, bin, sdpar);
424 }
425 }
426 }
427 }
428
429// delete ftmin;
430 for(lay=0; lay<MdcCalNLayer; lay++){
431 for(lr=0; lr<MdcCalLR; lr++){
432 delete ftmin[lay][lr];
433 delete ftmax[lay][lr];
434 }
435 }
436 return 1;
437}
void Fit()
Definition: Eangle1D/Fit.cxx:3
const int MdcCalNENTRSD
Definition: MdcCalParams.h:19
const int MdcCalSdNBIN
Definition: MdcCalParams.h:20
const int MdcCalNENTRXT
Definition: MdcCalParams.h:12
double initTm[MdcCalNLayer]
Definition: MdcCalParams.h:78
int fgCalib[MdcCalNLayer]
Definition: MdcCalParams.h:75
double t0Shift
Definition: MdcCalParams.h:50
double initT0
Definition: MdcCalParams.h:49
double tminFitRange[MdcCalNLayer][2]
Definition: MdcCalParams.h:76
double initSigma
Definition: MdcCalParams.h:53
double tminFitChindf
Definition: MdcCalParams.h:51
double tmaxFitChindf
Definition: MdcCalParams.h:52
double tmaxFitRange[MdcCalNLayer][2]
Definition: MdcCalParams.h:77
double getT0(int wireid) const
double getXtpar(int lay, int entr, int lr, int order)
int Id(void) const
Definition: MdcGeoWire.h:127
std::ofstream ofstream
Definition: bpkt_streams.h:42

The documentation for this class was generated from the following files: