BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
PreT0MdcCalib.cxx
Go to the documentation of this file.
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
14#include "EventModel/Event.h"
15#include "MdcRawEvent/MdcDigi.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
32 m_nzbin = 11;
33}
34
36}
37
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}
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
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
180 // dr cut
181 double dr = rectrk->getDr();
182 if(fabs(dr) > m_param.drCut) continue;
183
184 // dz cut
185 double dz = rectrk->getDz();
186 if(fabs(dz) > m_param.dzCut) continue;
187
188 // momentum cut
189 double p = rectrk -> getP();
190 if((fabs(p) < m_param.pCut[0]) || (fabs(p) > m_param.pCut[1])) continue;
191
192 // select cosmic track with y<0
193 double phi0 = rectrk -> getPhi0();
194 double phiTrk = phi0 + CLHEP::halfpi;
195 if(phiTrk > CLHEP::twopi) phiTrk -= CLHEP::twopi;
196 if(m_param.cosmicDwTrk && (phiTrk<CLHEP::pi)) continue;
197// cout << "cosmicDwTrk " << m_param.cosmicDwTrk << setw(15) << phiTrk*180./3.14 << endl;
198
199 nhit = rectrk -> getNHits();
200 for(k=0; k<nhit; k++){
201 rechit = rectrk -> getRecHit(k);
202 lay = rechit -> getLayid();
203 cel = rechit -> getCellid();
204 lr = rechit -> getLR();
205 zhit = rechit -> getZhit();
206 tdc = rechit -> getTdc();
207 tdr = rechit -> getTdrift();
208
209 traw = tdc * MdcCalTdcCnv;
210// traw = tdc; // ns
211
212 zprop = fabs(zhit - m_zst[lay]);
213 bin = (int)(zprop / m_zwid[lay]);
214 tp = zprop / m_vp[lay];
215 t0 = m_mdcFunSvc -> getT0(lay, cel);
216 tof = rechit -> getTof();
217
218 double adc = rechit -> getQhit();
219 double tw = m_mdcFunSvc->getTimeWalk(lay, adc);
220// cout << setw(15) << adc << setw(15) << tw << endl;
221
222 const MdcGeoWire* pWire = m_mdcGeomSvc -> Wire(lay, cel);
223 double y1 = pWire -> Backward().y();
224 double y2 = pWire -> Forward().y();
225 double ymid = (y1+y2)*0.5;
226// if(5==m_param.particle){ // cosmic-ray
227// if(ymid > 0) trec = traw - tp + tof - tes + m_param.timeShift - tw;
228// else trec = traw - tp - tof - tes + m_param.timeShift - tw;
229// } else{
230// trec = traw - tp - tof - tes + m_param.timeShift - tw;
231// }
232 trec = traw - tp - tof - tes + m_param.timeShift - tw;
233
234 m_hTrec[lay][lr] -> Fill(trec);
235 m_hTrec[lay][2] -> Fill(trec); // l-r in common
236 if(ymid > 0) m_hTrecCosm[lay][0] -> Fill(trec);
237 else m_hTrecCosm[lay][1] -> Fill(trec);
238
239 if(bin < m_nzbin){
240 m_hTrecZ[lay][lr][bin] -> Fill(trec);
241 m_hTrecZ[lay][2][bin] -> Fill(trec);
242 }
243 }
244 }
245 return 1;
246}
247
249}
250
252 IMessageSvc* msgSvc;
253 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
254 MsgStream log(msgSvc, "PreT0MdcCalib");
255 log << MSG::DEBUG << "PreT0MdcCalib::updateConst()" << endreq;
256
257// MdcCalib::updateConst(calconst);
258
259 // fit Tmin
260 int lay;
261 int wir;
262 int lr;
263 double t0Fit[MdcCalNLayer][MdcCalLR];
264 double t0Cal[MdcCalNLayer][MdcCalLR];
265 double tmax[MdcCalNLayer][MdcCalLR];
266 double initT0 = m_param.initT0;
267
268 int fitTminFg[MdcCalNLayer][MdcCalLR];
269 int fitTmaxFg[MdcCalNLayer][MdcCalLR];
270 double chisq;
271 double ndf;
272 double chindfTmin[MdcCalNLayer][MdcCalLR];
273 double chindfTmax[MdcCalNLayer][MdcCalLR];
274 char funname[200];
275
276 // add for cosmic-ray
277 TF1* ftminCosm[MdcCalNLayer][2];
278 double t0FitCosm[MdcCalNLayer][2];
279
280 TF1* ftmin[MdcCalNLayer][MdcCalLR];
281// sprintf(funname, "ftmin");
282// TF1* ftmin = new TF1(funname, funTmin, 0, 150, 6);
283 for(lay=0; lay<MdcCalNLayer; lay++){
284 for(lr=0; lr<MdcCalLR; lr++){
285 fitTminFg[lay][lr] = 0;
286 chindfTmin[lay][lr] = -1;
287 sprintf(funname, "ftmin%02d_%d", lay, lr);
288 ftmin[lay][lr] = new TF1(funname, funTmin, 0, 150, 6);
289
290 if(1 == m_param.fgCalib[lay]){
291 Stat_t nEntryTot = 0;
292 for(int ibin=1; ibin<=25; ibin++){
293 Stat_t entry = m_hTrec[lay][lr]->GetBinContent(ibin);
294 nEntryTot += entry;
295 }
296 double c0Ini = (double)nEntryTot / 25.0;
297 double c1Ini = (m_hTrec[lay][lr]->GetMaximum()) - c0Ini;
298
299 ftmin[lay][lr] -> SetParameter(0, c0Ini);
300 ftmin[lay][lr] -> SetParameter(1, c1Ini);
301 ftmin[lay][lr] -> SetParameter(2, 0);
302 ftmin[lay][lr] -> SetParameter(4, initT0);
303 ftmin[lay][lr] -> SetParameter(5, 2);
304
305 m_hTrec[lay][lr] -> Fit(funname, "Q", "",
306 m_param.tminFitRange[lay][0],
307 m_param.tminFitRange[lay][1]);
308 gStyle -> SetOptFit(11);
309// chisq = ftmin[lay][lr]->GetChisquare();
310// ndf = ftmin[lay][lr]->GetNDF();
311 chisq = ftmin[lay][lr]->GetChisquare();
312 ndf = ftmin[lay][lr]->GetNDF();
313 chindfTmin[lay][lr] = chisq / ndf;
314
315 if(chindfTmin[lay][lr] < m_param.tminFitChindf){
316 fitTminFg[lay][lr] = 1;
317 t0Fit[lay][lr] = ftmin[lay][lr]->GetParameter(4);
318// t0Fit[lay][lr] = ftmin->GetParameter(4);
319 t0Fit[lay][lr] += m_param.t0Shift;
320 t0Cal[lay][lr] = t0Fit[lay][lr] - m_param.timeShift;
321 }
322 }
323
324 if(0 == fitTminFg[lay][lr]){
325 wir = m_mdcGeomSvc->Wire(lay, 0)->Id();
326 t0Cal[lay][lr] = calconst->getT0(wir);
327 t0Fit[lay][lr] = t0Cal[lay][lr] + m_param.timeShift;
328 }
329 }
330
331 for(int iud=0; iud<2; iud++){
332 sprintf(funname, "ftminCosm_%02d_%d", lay, iud);
333 ftminCosm[lay][iud] = new TF1(funname, funTmin, 0, 150, 6);
334 ftminCosm[lay][iud] -> SetParameter(0, 0);
335 ftminCosm[lay][iud] -> SetParameter(4, initT0);
336 ftminCosm[lay][iud] -> SetParameter(5, 1);
337 m_hTrecCosm[lay][iud] -> Fit(funname, "Q", "",
338 m_param.tminFitRange[lay][0],
339 m_param.tminFitRange[lay][1]);
340 gStyle -> SetOptFit(11);
341 t0FitCosm[lay][iud] += m_param.t0Shift;
342 t0FitCosm[lay][iud] = ftminCosm[lay][iud]->GetParameter(4);
343 }
344 }
345
346 // fit Tmax
347 TF1* ftmax[MdcCalNLayer][MdcCalLR];
348 for(lay=0; lay<MdcCalNLayer; lay++){
349 for(lr=0; lr<MdcCalLR; lr++){
350 fitTmaxFg[lay][lr] = 0;
351 chindfTmax[lay][lr] = -1;
352 sprintf(funname, "ftmax%02d_%d", lay, lr);
353 ftmax[lay][lr] = new TF1(funname, funTmax, 250, 500, 4);
354
355 if(1 == m_param.fgCalib[lay]){
356 ftmax[lay][lr] -> SetParameter(2, m_param.initTm[lay]);
357 ftmax[lay][lr] -> SetParameter(3, 10);
358 m_hTrec[lay][lr] -> Fit(funname, "Q+", "",
359 m_param.tmaxFitRange[lay][0],
360 m_param.tmaxFitRange[lay][1]);
361 gStyle -> SetOptFit(11);
362 chisq = ftmax[lay][lr]->GetChisquare();
363 ndf = ftmax[lay][lr]->GetNDF();
364 chindfTmax[lay][lr] = chisq / ndf;
365 if(chindfTmax[lay][lr] < m_param.tmaxFitChindf){
366 fitTmaxFg[lay][lr] = 1;
367 tmax[lay][lr] = ftmax[lay][lr]->GetParameter(2);
368 }
369 }
370
371 if(0 == fitTmaxFg[lay][lr]){
372 tmax[lay][lr] = (calconst->getXtpar(lay, 0, lr, 6)) + t0Fit[lay][2];
373 }
374 }
375 }
376
377 // output for check
378 ofstream ft0("preT0.dat");
379 for(lay=0; lay<MdcCalNLayer; lay++){
380 ft0 << setw(5) << lay << setw(3) << fitTminFg[lay][2]
381 << setw(15) << t0Cal[lay][2] << setw(15) << t0Fit[lay][2]
382 << setw(15) << chindfTmin[lay][2] << endl;
383 }
384 ft0 << endl;
385 for(lay=0; lay<MdcCalNLayer; lay++){
386 ft0 << setw(5) << lay
387 << setw(3) << fitTmaxFg[lay][0] << setw(10) << tmax[lay][0]
388 << setw(10) << chindfTmax[lay][0]
389 << setw(3) << fitTmaxFg[lay][1] << setw(10) << tmax[lay][1]
390 << setw(10) << chindfTmax[lay][1]
391 << setw(3) << fitTmaxFg[lay][2] << setw(10) << tmax[lay][2]
392 << setw(10) << chindfTmax[lay][2]
393 << setw(10) << tmax[lay][0] - t0Fit[lay][2]
394 << setw(10) << tmax[lay][1] - t0Fit[lay][2]
395 << setw(10) << tmax[lay][2] - t0Fit[lay][2]
396 << endl;
397 }
398 ft0.close();
399 cout << "preT0.dat was written." << endl;
400
401 // output for cosmic T0
402 ofstream ft0cosm("cosmicT0.dat");
403 for(lay=0; lay<MdcCalNLayer; lay++){
404 ft0cosm << setw(5) << lay << setw(15) << t0Fit[lay][2]
405 << setw(15) << t0FitCosm[lay][0] << setw(15) << t0FitCosm[lay][1] << endl;
406 }
407 ft0cosm.close();
408
409 // set T0
410 int i;
411 int nwire = m_mdcGeomSvc -> getWireSize();
412 for(i=0; i<nwire; i++){
413 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
414 if(1 == m_param.fgCalib[lay]){
415 calconst -> resetT0(i, t0Cal[lay][2]);
416 calconst -> resetDelT0(i, 0.0);
417 }
418 }
419
420 // set tm of X-T
421 if(m_param.preT0SetTm){
422 int iEntr;
423 double tm;
424 for(lay=0; lay<MdcCalNLayer; lay++){
425 if(1 != m_param.fgCalib[lay]) continue;
426
427 for(iEntr=0; iEntr<MdcCalNENTRXT; iEntr++){
428 for(lr=0; lr<MdcCalLR; lr++){
429 tm = tmax[lay][lr] - t0Fit[lay][2];
430 if( (tmax[lay][lr] > m_param.tmaxFitRange[lay][0]) &&
431 (tmax[lay][lr] < m_param.tmaxFitRange[lay][1]) ){
432 calconst -> resetXtpar(lay, iEntr, lr, 6, tm);
433 }
434 }
435 }
436 }
437 }
438
439 // set sigma
440 int bin;
441 double sdpar = m_param.initSigma; // mm
442 for(lay=0; lay<MdcCalNLayer; lay++){
443 for(int iEntr=0; iEntr<MdcCalNENTRSD; iEntr++){
444 for(lr=0; lr<2; lr++){
445 for(bin=0; bin<MdcCalSdNBIN; bin++){
446 calconst -> resetSdpar(lay, iEntr, lr, bin, sdpar);
447 }
448 }
449 }
450 }
451
452// delete ftmin;
453 for(lay=0; lay<MdcCalNLayer; lay++){
454 for(lr=0; lr<MdcCalLR; lr++){
455 delete ftmin[lay][lr];
456 delete ftmax[lay][lr];
457 }
458 }
459 return 1;
460}
461Double_t PreT0MdcCalib::funTmin(Double_t* x, Double_t* par){
462 Double_t fitval;
463 fitval = par[0] + par[1]*exp( -par[2]*(x[0]-par[3]) ) /
464 ( 1 + exp( -(x[0]-par[4])/par[5] ));
465 return fitval;
466}
467
468Double_t PreT0MdcCalib::funTmax(Double_t* x, Double_t* par){
469 Double_t fitval;
470 fitval = par[0] + par[1] / (1 + exp((x[0]-par[2])/par[3]));
471 return fitval;
472}
473
curve Fill()
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)
void Fit()
Definition: Eangle1D/Fit.cxx:3
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
*******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 MdcCalNENTRSD
Definition: MdcCalParams.h:19
const double MdcCalTdcCnv
Definition: MdcCalParams.h:24
const int MdcCalSdNBIN
Definition: MdcCalParams.h:20
const int MdcCalNENTRXT
Definition: MdcCalParams.h:12
const int MdcCalLR
Definition: MdcCalParams.h:10
IMessageSvc * msgSvc()
#define NULL
virtual double getTimeWalk(int layid, double Q) const =0
virtual const MdcGeoWire *const Wire(unsigned id)=0
bool cosmicDwTrk
Definition: MdcCalParams.h:71
double initTm[MdcCalNLayer]
Definition: MdcCalParams.h:83
int fgCalib[MdcCalNLayer]
Definition: MdcCalParams.h:80
double t0Shift
Definition: MdcCalParams.h:51
double pCut[2]
Definition: MdcCalParams.h:75
double initT0
Definition: MdcCalParams.h:50
double tminFitRange[MdcCalNLayer][2]
Definition: MdcCalParams.h:81
double initSigma
Definition: MdcCalParams.h:54
double dzCut
Definition: MdcCalParams.h:74
double tminFitChindf
Definition: MdcCalParams.h:52
double timeShift
Definition: MdcCalParams.h:43
double tmaxFitChindf
Definition: MdcCalParams.h:53
double drCut
Definition: MdcCalParams.h:73
double tmaxFitRange[MdcCalNLayer][2]
Definition: MdcCalParams.h:82
double getDz() const
Definition: MdcCalRecTrk.h:33
double getDr() const
Definition: MdcCalRecTrk.h:30
double getT0(int wireid) const
double getXtpar(int lay, int entr, int lr, int order)
virtual void clear()=0
Definition: MdcCalib.cxx:80
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
int Id(void) const
Definition: MdcGeoWire.h:127
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128
int updateConst(MdcCalibConst *calconst)
int fillHist(MdcCalEvent *event)
void initialize(TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)
void printCut() const