BOSS 7.0.3
BESIII Offline Software System
Loading...
Searching...
No Matches
Wr2dMdcCalib.cxx
Go to the documentation of this file.
1#include "MdcCalibAlg/Wr2dMdcCalib.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 "TF1.h"
10#include "TMinuit.h"
11
12#include <iostream>
13#include <fstream>
14#include <iomanip>
15#include <cstring>
16
17using namespace std;
18
19bool Wr2dMdcCalib::fgBIN[MdcCalWrNBin];
20double Wr2dMdcCalib::xBIN[MdcCalWrNBin];
21double Wr2dMdcCalib::yBIN[MdcCalWrNBin];
22double Wr2dMdcCalib::zBIN[MdcCalWrNBin];
23double Wr2dMdcCalib::zBINERR[MdcCalWrNBin];
24double Wr2dMdcCalib::zMIN;
25double Wr2dMdcCalib::zMAX;
26
27Wr2dMdcCalib::Wr2dMdcCalib(){
28}
29
30Wr2dMdcCalib::~Wr2dMdcCalib(){
31}
32
33void Wr2dMdcCalib::clear(){
34 int bin;
35 for(int i=0; i<MdcCalTotCell; i++){
36 for(bin=0; bin<MdcCalWrNBin; bin++){
37 delete m_hl[i][bin];
38 delete m_hr[i][bin];
39 }
40 }
41 delete m_fdWire;
42
43 MdcCalib::clear();
44}
45
46void Wr2dMdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
47 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
48 IMessageSvc *msgSvc;
49 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
50 MsgStream log(msgSvc, "Wr2dMdcCalib");
51 log << MSG::INFO << "Wr2dMdcCalib::initialize()" << endreq;
52
53 m_hlist = hlist;
54 m_mdcGeomSvc = mdcGeomSvc;
55 m_mdcFunSvc = mdcFunSvc;
56 m_mdcUtilitySvc = mdcUtilitySvc;
57
58 MdcCalib::initialize(m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc);
59
60 int i;
61 int bin;
62 int lay;
63 int cel;
64 char hname[200];
65
66 m_fdWire = new TFolder("WireCor", "WireCor");
67 m_hlist->Add(m_fdWire);
68
69 for(i=0; i<MdcCalTotCell; i++){
70 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
71 cel = m_mdcGeomSvc -> Wire(i) -> Cell();
72
73 for(bin=0; bin<MdcCalWrNBin; bin++){
74 sprintf(hname, "h%04d_L_%02d", i, bin);
75 m_hl[i][bin] = new TH1F(hname, "", 300, -1.5, 1.5);
76 m_fdWire->Add(m_hl[i][bin]);
77
78 sprintf(hname, "h%04d_R_%02d", i, bin);
79 m_hr[i][bin] = new TH1F(hname, "", 300, -1.5, 1.5);
80 m_fdWire->Add(m_hr[i][bin]);
81 }
82 }
83
84 for(lay=0; lay<MdcCalNLayer; lay++){
85 m_zwest[lay] = m_mdcGeomSvc->Wire(lay, 0)->Forward().z();
86 m_zeast[lay] = m_mdcGeomSvc->Wire(lay, 0)->Backward().z();
87 m_zwid[lay] = (m_zeast[lay] - m_zwest[lay]) / (double)MdcCalWrNBin;
88
89 for(bin=0; bin<MdcCalWrNBin; bin++){
90 m_zbinCen[lay][bin] = ((double)bin + 0.5) * m_zwid[lay] + m_zwest[lay];
91 }
92 }
93}
94
95int Wr2dMdcCalib::fillHist(MdcCalEvent* event){
96 IMessageSvc *msgSvc;
97 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
98 MsgStream log(msgSvc, "Wr2dMdcCalib");
99 log << MSG::DEBUG << "Wr2dMdcCalib::fillHist()" << endreq;
100
101 MdcCalib::fillHist(event);
102
103 // get EsTimeCol
104 bool esCutFg = event->getEsCutFlag();
105 if( ! esCutFg ) return -1;
106
107 int i;
108 int k;
109 int ntrk;
110 int nhit;
111
112 int bin;
113 int lay;
114 int cel;
115 int wir;
116 int lr;
117 double dmeas;
118 double doca;
119 double residual;
120 double zhit;
121
122 MdcCalRecTrk* rectrk;
123 MdcCalRecHit* rechit;
124
125 ntrk = event->getNTrk();
126 log << MSG::DEBUG << "number of tracks: " << ntrk << endreq;
127
128 for(i=0; i<ntrk; i++){
129 rectrk = event -> getRecTrk(i);
130 nhit = rectrk -> getNHits();
131
132 log << MSG::DEBUG << "number of hits: " << nhit << endreq;
133 for(k=0; k<nhit; k++){
134 rechit = rectrk -> getRecHit(k);
135 lay = rechit -> getLayid();
136 cel = rechit -> getCellid();
137 wir = m_mdcGeomSvc -> Wire(lay, cel) -> Id();
138 lr = rechit -> getLR();
139 doca = rechit -> getDocaInc();
140 dmeas = rechit -> getDmeas();
141 residual = rechit -> getResiInc();
142 zhit = rechit -> getZhit();
143
144 if((wir<MdcCalTotCell) && (fabs(zhit) < m_zeast[lay])){
145 bin = (int)(zhit / m_zwid[lay]);
146 if( 0 == lr ){
147 m_hl[wir][bin] -> Fill(residual);
148 }else if( 1 == lr ){
149 m_hr[wir][bin] -> Fill(residual);
150 }
151 }else{
152 std::cout << "wir: " << wir << std::endl;
153 }
154 }
155 }
156 return 1;
157}
158
159int Wr2dMdcCalib::updateConst(MdcCalibConst* calconst){
160 IMessageSvc *msgSvc;
161 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
162 MsgStream log(msgSvc, "Wr2dMdcCalib");
163 log << MSG::INFO << "Wr2dMdcCalib::updateConst()" << endreq;
164
165 int i;
166 int k;
167 int bin;
168 int lay;
169 int cel;
170 double dw;
171
172 // minuit for wires
173 Int_t ierflg;
174 Int_t istat;
175 Int_t nvpar;
176 Int_t nparx;
177 Double_t fmin;
178 Double_t edm;
179 Double_t errdef;
180 Double_t arglist[10];
181
182 TMinuit *gmtrw = new TMinuit(5);
183 gmtrw -> SetPrintLevel(-1);
184 gmtrw -> SetFCN(fcnWireParab);
185 gmtrw -> SetErrorDef(1.0);
186 gmtrw -> mnparm(0, "xf", 0.0, 0.01, 0, 0, ierflg);
187 gmtrw -> mnparm(1, "yf", 0.0, 0.01, 0, 0, ierflg);
188 gmtrw -> mnparm(2, "xb", 0.0, 0.01, 0, 0, ierflg);
189 gmtrw -> mnparm(3, "yb", 0.0, 0.01, 0, 0, ierflg);
190 gmtrw -> mnparm(4, "ten", 0.0, 0.1, 0, 0, ierflg);
191 arglist[0] = 0;
192 gmtrw -> mnexcm("Set NOW", arglist, 0, ierflg);
193
194 double a;
195 double dx;
196 double dy;
197 double dz;
198 double length;
199 double ten[] = {15, 15, 15, 16, 16, 17, 17, 18, 14, 14,
200 19, 19, 24, 24, 31, 31, 37, 37, 45, 45,
201 46, 47, 47, 47, 47, 48, 48, 48, 48, 49,
202 49, 49, 49, 50, 50, 50, 51, 51, 51, 52,
203 52, 52, 52};
204 double wpar[5];
205 double wparErr[5];
206 double sinphiw;
207 double cosphiw;
208 double deldw;
209 double delxw;
210 double delyw;
211
212 int nFit;
213 Stat_t entryL;
214 Stat_t entryR;
215 Double_t hcen;
216 Double_t cenL[MdcCalWrNBin];
217 Double_t errL[MdcCalWrNBin];
218 Double_t cenR[MdcCalWrNBin];
219 Double_t errR[MdcCalWrNBin];
220
221 ofstream fwlog("logWireCor.txt");
222 ofstream fwpc("wireposCor.txt");
223 ofstream fwpcErr("wireposErr.txt");
224
225 fwpc << setw(5) << "wireId" << setw(15) << "dx_East(mm)"
226 << setw(15) << "dy_East(mm)" << setw(15) << "dz_East(mm)"
227 << setw(15) << "dx_West(mm)" << setw(15) << "dy_West(mm)"
228 << setw(15) << "dz_West(mm)" << endl;
229
230 for(i=0; i<MdcCalTotCell; i++){
231 for(k=0; k<5; k++) wparErr[k] = 999.0;
232 nFit = 0;
233 for(bin=0; bin<MdcCalWrNBin; bin++){
234 entryL = m_hl[i][bin] -> GetEntries();
235 entryR = m_hr[i][bin] -> GetEntries();
236 if((entryL < 100) && (entryR < 100)){
237 fgBIN[bin] = false;
238 continue;
239 } else{
240 fgBIN[bin] = true;
241 }
242 nFit++;
243
244 if(1 == m_param.fgCalib[lay]){
245 hcen = m_hl[i][bin] -> GetMean();
246 if(entryL > 500){
247 m_hl[i][bin] -> Fit("gaus", "Q", "", hcen-1.5, hcen+1.5);
248 cenL[bin] = m_hl[i][bin]->GetFunction("gaus")->GetParameter(1);
249 errL[bin] = m_hl[i][bin]->GetFunction("gaus")->GetParameter(2);
250 }
251 else{
252 cenL[bin] = hcen;
253 errL[bin] = m_hl[i][bin] -> GetRMS();
254 }
255
256 hcen = m_hr[i][bin] -> GetMean();
257 if(entryR > 500){
258 m_hr[i][bin] -> Fit("gaus", "Q", "", hcen-1.5, hcen+1.5);
259 cenR[bin] = m_hr[i][bin]->GetFunction("gaus")->GetParameter(1);
260 errR[bin] = m_hr[i][bin]->GetFunction("gaus")->GetParameter(2);
261 }
262 else{
263 cenR[bin] = hcen;
264 errR[bin] = m_hr[i][bin] -> GetRMS();
265 }
266 } else{
267 cenL[bin] = 0.0;
268 errL[bin] = 0.15;
269 cenR[bin] = 0.0;
270 errR[bin] = 0.15;
271 }
272 }
273 if(nFit < 8) continue;
274
275 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
276 cel = m_mdcGeomSvc -> Wire(i) -> Cell();
277 zMIN = m_zwest[lay];
278 zMAX = m_zeast[lay];
279 wpar[0] = m_mdcGeomSvc->Wire(lay, 0)->Backward().x();
280 wpar[1] = m_mdcGeomSvc->Wire(lay, 0)->Backward().y();
281 wpar[2] = m_mdcGeomSvc->Wire(lay, 0)->Forward().x();
282 wpar[3] = m_mdcGeomSvc->Wire(lay, 0)->Forward().y();
283 wpar[4] = ten[lay];
284
285 a = 9.47e-06 / (2 * wpar[4]);
286 dx = wpar[0] - wpar[2];
287 dy = wpar[1] - wpar[3];
288 dz = zMAX - zMIN;
289 length = sqrt(dx*dx + dz*dz);
290
291 for(bin=0; bin<MdcCalWrNBin; bin++){
292 zBIN[bin] = m_zbinCen[lay][bin];
293 xBIN[bin] = (wpar[0]-wpar[2])*(zBIN[bin]-zMIN)/(zMAX-zMIN) + wpar[2];
294 yBIN[bin] = a*zBIN[bin]*zBIN[bin] + (wpar[1]-wpar[3])*zBIN[bin]/length
295 + 0.5*(wpar[1]+wpar[3]) - a*length*length/4.0;
296
297 sinphiw = yBIN[bin] / sqrt(xBIN[bin]*xBIN[bin] + yBIN[bin]*yBIN[bin]);
298 cosphiw = xBIN[bin] / sqrt(xBIN[bin]*xBIN[bin] + yBIN[bin]*yBIN[bin]);
299
300 deldw = - (cenL[bin]-cenR[bin])/2.0;
301 delxw = -deldw * sinphiw;
302 delyw = deldw * cosphiw;
303
304 fwlog << setw(3) << lay << setw(4) << cel << setw(5) << i
305 << setw(4) << bin << setw(15) << zBIN[bin]
306 << setw(15) << deldw << setw(15) << delxw
307 << setw(15) << delyw << endl;
308
309 xBIN[bin] += delxw;
310 yBIN[bin] += delyw;
311
312 zBINERR[bin] = sqrt((errL[bin]*errL[bin]) + (errR[bin]*errR[bin]))/2.0;
313 }
314
315 arglist[0] = 1;
316 arglist[1] = wpar[0];
317 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
318 arglist[0] = 2;
319 arglist[1] = wpar[1];
320 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
321 arglist[0] = 3;
322 arglist[1] = wpar[2];
323 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
324 arglist[0] = 4;
325 arglist[1] = wpar[3];
326 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
327 arglist[0] = 5;
328 arglist[1] = wpar[4];
329 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
330 gmtrw -> mnexcm("FIX", arglist, 1, ierflg);
331
332 arglist[0] = 2000;
333 arglist[1] = 0.1;
334 gmtrw -> mnexcm("MIGRAD", arglist, 2, ierflg);
335 gmtrw -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
336
337 if( (0==ierflg) && (3==istat) ){
338 gmtrw -> GetParameter(0, wpar[0], wparErr[0]);
339 gmtrw -> GetParameter(1, wpar[1], wparErr[1]);
340 gmtrw -> GetParameter(2, wpar[2], wparErr[2]);
341 gmtrw -> GetParameter(3, wpar[3], wparErr[3]);
342 gmtrw -> GetParameter(4, wpar[4], wparErr[4]);
343 }
344 gmtrw -> mnexcm("RELease", arglist, 1, ierflg);
345
346 fwlog << setw(5) << i << setw(15) << wpar[0] << setw(15) << wpar[1]
347 << setw(15) << wpar[2] << setw(15) << wpar[3] << endl;
348 fwpc << setw(5) << i << setw(15) << wpar[0] << setw(15) << wpar[1]
349 << setw(15) << "0" << setw(15) << wpar[2]
350 << setw(15) << wpar[3] << setw(15) << "0" << endl;
351 fwpcErr << setw(5) << i << setw(15) << wparErr[0] << setw(15) << wparErr[1]
352 << setw(15) << wparErr[2] << setw(15) << wparErr[3] << endl;
353 }
354 fwlog.close();
355 fwpc.close();
356 fwpcErr.close();
357
358 delete gmtrw;
359 return 1;
360}
361
362void Wr2dMdcCalib::fcnWireParab(Int_t &npar, Double_t *gin,
363 Double_t &f, Double_t *par, Int_t iflag){
364 Int_t bin;
365 Double_t xfit;
366 Double_t yfit;
367
368 double a = 9.47e-06 / (2 * par[4]);
369 double dx = par[0] - par[2];
370 double dy = par[1] - par[3];
371 double dz = zMAX - zMIN;
372 double length = sqrt(dx*dx + dz*dz);
373
374 Double_t chisq = 0.0;
375 Double_t deta;
376
377 for(bin=0; bin<MdcCalWrNBin; bin++){
378 if( fgBIN[bin] ){
379 xfit = (par[0] - par[2]) * (zBIN[bin] - zMIN) / (zMAX - zMIN) + par[2];
380 yfit = a*zBIN[bin]*zBIN[bin] + (par[1] - par[3])*zBIN[bin]/length
381 + 0.5*(par[1] + par[3]) - a*length*length/4.0;
382
383 deta = ( (xfit-xBIN[bin])*(xfit-xBIN[bin])+
384 (yfit-yBIN[bin])*(yfit-yBIN[bin]) ) / (zBINERR[bin]*zBINERR[bin]);
385 chisq += deta;
386 }
387 }
388 f = chisq;
389}
390
Fit
gr Fit("g1","R")
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
IMdcGeomSvc::Wire
virtual const MdcGeoWire *const Wire(unsigned id)=0
MdcCalib::clear
virtual void clear()=0
Definition: MdcCalib.cxx:76
MdcCalib::initialize
virtual void initialize(TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)=0
Definition: MdcCalib.cxx:202
MdcCalib::fillHist
virtual int fillHist(MdcCalEvent *event)=0
Definition: MdcCalib.cxx:692
Wr2dMdcCalib::fillHist
int fillHist(MdcCalEvent *event)
Definition: Wr2dMdcCalib.cxx:95
Wr2dMdcCalib::fcnWireParab
static void fcnWireParab(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: Wr2dMdcCalib.cxx:362
Wr2dMdcCalib::initialize
void initialize(TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)
Definition: Wr2dMdcCalib.cxx:46
Wr2dMdcCalib::updateConst
int updateConst(MdcCalibConst *calconst)
Definition: Wr2dMdcCalib.cxx:159