BOSS 6.6.4.p01
BESIII Offline Software System
Loading...
Searching...
No Matches
PreT0MdcCalib.cxx
Go to the documentation of this file.
1#include "MdcCalibAlg/PreT0MdcCalib.h"
2
3#include "GaudiKernel/MsgStream.h"
4#include "GaudiKernel/IMessageSvc.h"
5#include "GaudiKernel/StatusCode.h"
6#include "GaudiKernel/ISvcLocator.h"
7#include "GaudiKernel/Bootstrap.h"
8
9#include "GaudiKernel/SmartDataPtr.h"
10#include "GaudiKernel/IDataProviderSvc.h"
11#include "GaudiKernel/PropertyMgr.h"
12
13#include "EvTimeEvent/RecEsTime.h"
14#include "EventModel/Event.h"
15#include "MdcRawEvent/MdcDigi.h"
16#include "Identifier/Identifier.h"
17#include "Identifier/MdcID.h"
18#include "TStyle.h"
19
20#include <iostream>
21#include <fstream>
22#include <iomanip>
23#include <cstring>
24
25#include <math.h>
26#include "TH2F.h"
27#include "TH1F.h"
28#include "TFolder.h"
29#include "TF1.h"
30
31PreT0MdcCalib::PreT0MdcCalib(){
32 m_nzbin = 11;
33}
34
35PreT0MdcCalib::~PreT0MdcCalib(){
36}
37
38void PreT0MdcCalib::clear(){
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
56 MdcCalib::clear();
57}
58
59void PreT0MdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
60 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
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}
137
138int PreT0MdcCalib::fillHist(MdcCalEvent* event){
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}
227
228int PreT0MdcCalib::updateConst(MdcCalibConst* calconst){
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}
438Double_t PreT0MdcCalib::funTmin(Double_t* x, Double_t* par){
439 Double_t fitval;
440 fitval = par[0] + par[1]*exp( -par[2]*(x[0]-par[3]) ) /
441 ( 1 + exp( -(x[0]-par[4])/par[5] ));
442 return fitval;
443}
444
445Double_t PreT0MdcCalib::funTmax(Double_t* x, Double_t* par){
446 Double_t fitval;
447 fitval = par[0] + par[1] / (1 + exp((x[0]-par[2])/par[3]));
448 return fitval;
449}
450
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
MdcCalNLayer
const int MdcCalNLayer
Definition: MdcCalParams.h:6
MdcCalNENTRSD
const int MdcCalNENTRSD
Definition: MdcCalParams.h:19
MdcCalTdcCnv
const double MdcCalTdcCnv
Definition: MdcCalParams.h:24
MdcCalSdNBIN
const int MdcCalSdNBIN
Definition: MdcCalParams.h:20
MdcCalNENTRXT
const int MdcCalNENTRXT
Definition: MdcCalParams.h:12
MdcCalLR
const int MdcCalLR
Definition: MdcCalParams.h:10
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
IMdcCalibFunSvc::getTimeWalk
virtual double getTimeWalk(int layid, double Q) const =0
IMdcGeomSvc::Wire
virtual const MdcGeoWire *const Wire(unsigned id)=0
MdcCalParams::initTm
double initTm[MdcCalNLayer]
Definition: MdcCalParams.h:78
MdcCalParams::fgCalib
int fgCalib[MdcCalNLayer]
Definition: MdcCalParams.h:75
MdcCalParams::t0Shift
double t0Shift
Definition: MdcCalParams.h:50
MdcCalParams::initT0
double initT0
Definition: MdcCalParams.h:49
MdcCalParams::tminFitRange
double tminFitRange[MdcCalNLayer][2]
Definition: MdcCalParams.h:76
MdcCalParams::initSigma
double initSigma
Definition: MdcCalParams.h:53
MdcCalParams::tminFitChindf
double tminFitChindf
Definition: MdcCalParams.h:51
MdcCalParams::timeShift
double timeShift
Definition: MdcCalParams.h:42
MdcCalParams::tmaxFitChindf
double tmaxFitChindf
Definition: MdcCalParams.h:52
MdcCalParams::tmaxFitRange
double tmaxFitRange[MdcCalNLayer][2]
Definition: MdcCalParams.h:77
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
MdcCalib::clear
virtual void clear()=0
Definition: MdcCalib.cxx:76
MdcGeoWire::Forward
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
MdcGeoWire::Id
int Id(void) const
Definition: MdcGeoWire.h:127
MdcGeoWire::Backward
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128
PreT0MdcCalib::updateConst
int updateConst(MdcCalibConst *calconst)
Definition: PreT0MdcCalib.cxx:228
PreT0MdcCalib::fillHist
int fillHist(MdcCalEvent *event)
Definition: PreT0MdcCalib.cxx:138
PreT0MdcCalib::initialize
void initialize(TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)
Definition: PreT0MdcCalib.cxx:59