163 {
165 Gaudi::svcLocator() -> service(
"MessageSvc",
msgSvc);
166 MsgStream log(
msgSvc,
"Wr2dMdcCalib");
167 log << MSG::INFO << "Wr2dMdcCalib::updateConst()" << endreq;
168
169 int i;
170 int k;
172 int lay;
173 int cel;
174 double dw;
175
176
177 Int_t ierflg;
178 Int_t istat;
179 Int_t nvpar;
180 Int_t nparx;
181 Double_t fmin;
182 Double_t edm;
183 Double_t errdef;
184 Double_t arglist[10];
185
186 TMinuit *gmtrw = new TMinuit(5);
187 gmtrw -> SetPrintLevel(-1);
189 gmtrw -> SetErrorDef(1.0);
190 gmtrw -> mnparm(0, "xf", 0.0, 0.01, 0, 0, ierflg);
191 gmtrw -> mnparm(1, "yf", 0.0, 0.01, 0, 0, ierflg);
192 gmtrw -> mnparm(2, "xb", 0.0, 0.01, 0, 0, ierflg);
193 gmtrw -> mnparm(3, "yb", 0.0, 0.01, 0, 0, ierflg);
194 gmtrw -> mnparm(4, "ten", 0.0, 0.1, 0, 0, ierflg);
195 arglist[0] = 0;
196 gmtrw -> mnexcm("Set NOW", arglist, 0, ierflg);
197
198 double a;
199 double dx;
200 double dy;
201 double dz;
202 double length;
203 double ten[] = {15, 15, 15, 16, 16, 17, 17, 18, 14, 14,
204 19, 19, 24, 24, 31, 31, 37, 37, 45, 45,
205 46, 47, 47, 47, 47, 48, 48, 48, 48, 49,
206 49, 49, 49, 50, 50, 50, 51, 51, 51, 52,
207 52, 52, 52};
208 double wpar[5];
209 double wparErr[5];
210 double sinphiw;
211 double cosphiw;
212 double deldw;
213 double delxw;
214 double delyw;
215
216 int nFit;
217 Stat_t entryL;
218 Stat_t entryR;
219 Double_t hcen;
224
225 ofstream fwlog("logWireCor.txt");
226 ofstream fwpc("wireposCor.txt");
227 ofstream fwpcErr("wireposErr.txt");
228
229 fwpc << setw(5) << "wireId" << setw(15) << "dx_East(mm)"
230 << setw(15) << "dy_East(mm)" << setw(15) << "dz_East(mm)"
231 << setw(15) << "dx_West(mm)" << setw(15) << "dy_West(mm)"
232 << setw(15) << "dz_West(mm)" << endl;
233
235 for(k=0; k<5; k++) wparErr[k] = 999.0;
236 nFit = 0;
238 entryL = m_hl[i][
bin] -> GetEntries();
239 entryR = m_hr[i][
bin] -> GetEntries();
240 if((entryL < 100) && (entryR < 100)){
242 continue;
243 } else{
245 }
246 nFit++;
247
249 hcen = m_hl[i][
bin] -> GetMean();
250 if(entryL > 500){
251 m_hl[i][
bin] ->
Fit(
"gaus",
"Q",
"", hcen-1.5, hcen+1.5);
252 cenL[
bin] = m_hl[i][
bin]->GetFunction(
"gaus")->GetParameter(1);
253 errL[
bin] = m_hl[i][
bin]->GetFunction(
"gaus")->GetParameter(2);
254 }
255 else{
257 errL[
bin] = m_hl[i][
bin] -> GetRMS();
258 }
259
260 hcen = m_hr[i][
bin] -> GetMean();
261 if(entryR > 500){
262 m_hr[i][
bin] ->
Fit(
"gaus",
"Q",
"", hcen-1.5, hcen+1.5);
263 cenR[
bin] = m_hr[i][
bin]->GetFunction(
"gaus")->GetParameter(1);
264 errR[
bin] = m_hr[i][
bin]->GetFunction(
"gaus")->GetParameter(2);
265 }
266 else{
268 errR[
bin] = m_hr[i][
bin] -> GetRMS();
269 }
270 } else{
275 }
276 }
277 if(nFit < 8) continue;
278
279 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
280 cel = m_mdcGeomSvc -> Wire(i) -> Cell();
285 wpar[2] = m_mdcGeomSvc->
Wire(lay, 0)->
Forward().x();
286 wpar[3] = m_mdcGeomSvc->
Wire(lay, 0)->
Forward().y();
287 wpar[4] = ten[lay];
288
289 a = 9.47e-06 / (2 * wpar[4]);
290 dx = wpar[0] - wpar[2];
291 dy = wpar[1] - wpar[3];
293 length = sqrt(dx*dx + dz*dz);
294
299 + 0.5*(wpar[1]+wpar[3]) - a*length*length/4.0;
300
303
304 deldw = - (cenL[
bin]-cenR[
bin])/2.0;
305 delxw = -deldw * sinphiw;
306 delyw = deldw * cosphiw;
307
308 fwlog << setw(3) << lay << setw(4) << cel << setw(5) << i
310 << setw(15) << deldw << setw(15) << delxw
311 << setw(15) << delyw << endl;
312
315
317 }
318
319 arglist[0] = 1;
320 arglist[1] = wpar[0];
321 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
322 arglist[0] = 2;
323 arglist[1] = wpar[1];
324 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
325 arglist[0] = 3;
326 arglist[1] = wpar[2];
327 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
328 arglist[0] = 4;
329 arglist[1] = wpar[3];
330 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
331 arglist[0] = 5;
332 arglist[1] = wpar[4];
333 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
334 gmtrw -> mnexcm("FIX", arglist, 1, ierflg);
335
336 arglist[0] = 2000;
337 arglist[1] = 0.1;
338 gmtrw -> mnexcm("MIGRAD", arglist, 2, ierflg);
339 gmtrw -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
340
341 if( (0==ierflg) && (3==istat) ){
342 gmtrw -> GetParameter(0, wpar[0], wparErr[0]);
343 gmtrw -> GetParameter(1, wpar[1], wparErr[1]);
344 gmtrw -> GetParameter(2, wpar[2], wparErr[2]);
345 gmtrw -> GetParameter(3, wpar[3], wparErr[3]);
346 gmtrw -> GetParameter(4, wpar[4], wparErr[4]);
347 }
348 gmtrw -> mnexcm("RELease", arglist, 1, ierflg);
349
350 fwlog << setw(5) << i << setw(15) << wpar[0] << setw(15) << wpar[1]
351 << setw(15) << wpar[2] << setw(15) << wpar[3] << endl;
352 fwpc << setw(5) << i << setw(15) << wpar[0] << setw(15) << wpar[1]
353 << setw(15) << "0" << setw(15) << wpar[2]
354 << setw(15) << wpar[3] << setw(15) << "0" << endl;
355 fwpcErr << setw(5) << i << setw(15) << wparErr[0] << setw(15) << wparErr[1]
356 << setw(15) << wparErr[2] << setw(15) << wparErr[3] << endl;
357 }
358 fwlog.close();
359 fwpc.close();
360 fwpcErr.close();
361
362 delete gmtrw;
363 return 1;
364}
int fgCalib[MdcCalNLayer]
static void fcnWireParab(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)