159 {
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;
168 int lay;
169 int cel;
170 double dw;
171
172
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);
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;
220
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
231 for(k=0; k<5; k++) wparErr[k] = 999.0;
232 nFit = 0;
234 entryL = m_hl[i][
bin] -> GetEntries();
235 entryR = m_hr[i][
bin] -> GetEntries();
236 if((entryL < 100) && (entryR < 100)){
238 continue;
239 } else{
241 }
242 nFit++;
243
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{
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{
264 errR[
bin] = m_hr[i][
bin] -> GetRMS();
265 }
266 } else{
271 }
272 }
273 if(nFit < 8) continue;
274
275 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
276 cel = m_mdcGeomSvc -> Wire(i) -> Cell();
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];
289 length = sqrt(dx*dx + dz*dz);
290
295 + 0.5*(wpar[1]+wpar[3]) - a*length*length/4.0;
296
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
306 << setw(15) << deldw << setw(15) << delxw
307 << setw(15) << delyw << endl;
308
311
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}
int fgCalib[MdcCalNLayer]
static void fcnWireParab(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)