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