00001 #include "MdcCalibAlg/Wr2dMdcCalib.h"
00002
00003 #include "GaudiKernel/MsgStream.h"
00004 #include "GaudiKernel/IMessageSvc.h"
00005 #include "GaudiKernel/StatusCode.h"
00006 #include "GaudiKernel/ISvcLocator.h"
00007 #include "GaudiKernel/Bootstrap.h"
00008
00009 #include "TF1.h"
00010 #include "TMinuit.h"
00011
00012 #include <iostream>
00013 #include <fstream>
00014 #include <iomanip>
00015 #include <cstring>
00016
00017 using namespace std;
00018
00019 bool Wr2dMdcCalib::fgBIN[MdcCalWrNBin];
00020 double Wr2dMdcCalib::xBIN[MdcCalWrNBin];
00021 double Wr2dMdcCalib::yBIN[MdcCalWrNBin];
00022 double Wr2dMdcCalib::zBIN[MdcCalWrNBin];
00023 double Wr2dMdcCalib::zBINERR[MdcCalWrNBin];
00024 double Wr2dMdcCalib::zMIN;
00025 double Wr2dMdcCalib::zMAX;
00026
00027 Wr2dMdcCalib::Wr2dMdcCalib(){
00028 }
00029
00030 Wr2dMdcCalib::~Wr2dMdcCalib(){
00031 }
00032
00033 void Wr2dMdcCalib::clear(){
00034 int bin;
00035 for(int i=0; i<MdcCalTotCell; i++){
00036 for(bin=0; bin<MdcCalWrNBin; bin++){
00037 delete m_hl[i][bin];
00038 delete m_hr[i][bin];
00039 }
00040 }
00041 delete m_fdWire;
00042
00043 MdcCalib::clear();
00044 }
00045
00046 void Wr2dMdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
00047 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
00048 IMessageSvc *msgSvc;
00049 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
00050 MsgStream log(msgSvc, "Wr2dMdcCalib");
00051 log << MSG::INFO << "Wr2dMdcCalib::initialize()" << endreq;
00052
00053 m_hlist = hlist;
00054 m_mdcGeomSvc = mdcGeomSvc;
00055 m_mdcFunSvc = mdcFunSvc;
00056 m_mdcUtilitySvc = mdcUtilitySvc;
00057
00058 MdcCalib::initialize(m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc);
00059
00060 int i;
00061 int bin;
00062 int lay;
00063 int cel;
00064 char hname[200];
00065
00066 m_fdWire = new TFolder("WireCor", "WireCor");
00067 m_hlist->Add(m_fdWire);
00068
00069 for(i=0; i<MdcCalTotCell; i++){
00070 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
00071 cel = m_mdcGeomSvc -> Wire(i) -> Cell();
00072
00073 for(bin=0; bin<MdcCalWrNBin; bin++){
00074 sprintf(hname, "h%04d_L_%02d", i, bin);
00075 m_hl[i][bin] = new TH1F(hname, "", 300, -1.5, 1.5);
00076 m_fdWire->Add(m_hl[i][bin]);
00077
00078 sprintf(hname, "h%04d_R_%02d", i, bin);
00079 m_hr[i][bin] = new TH1F(hname, "", 300, -1.5, 1.5);
00080 m_fdWire->Add(m_hr[i][bin]);
00081 }
00082 }
00083
00084 for(lay=0; lay<MdcCalNLayer; lay++){
00085 m_zwest[lay] = m_mdcGeomSvc->Wire(lay, 0)->Forward().z();
00086 m_zeast[lay] = m_mdcGeomSvc->Wire(lay, 0)->Backward().z();
00087 m_zwid[lay] = (m_zeast[lay] - m_zwest[lay]) / (double)MdcCalWrNBin;
00088
00089 for(bin=0; bin<MdcCalWrNBin; bin++){
00090 m_zbinCen[lay][bin] = ((double)bin + 0.5) * m_zwid[lay] + m_zwest[lay];
00091 }
00092 }
00093 }
00094
00095 int Wr2dMdcCalib::fillHist(MdcCalEvent* event){
00096 IMessageSvc *msgSvc;
00097 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00098 MsgStream log(msgSvc, "Wr2dMdcCalib");
00099 log << MSG::DEBUG << "Wr2dMdcCalib::fillHist()" << endreq;
00100
00101 MdcCalib::fillHist(event);
00102
00103
00104 bool esCutFg = event->getEsCutFlag();
00105 if( ! esCutFg ) return -1;
00106
00107 int i;
00108 int k;
00109 int ntrk;
00110 int nhit;
00111
00112 int bin;
00113 int lay;
00114 int cel;
00115 int wir;
00116 int lr;
00117 double dmeas;
00118 double doca;
00119 double residual;
00120 double zhit;
00121
00122 MdcCalRecTrk* rectrk;
00123 MdcCalRecHit* rechit;
00124
00125 ntrk = event->getNTrk();
00126 log << MSG::DEBUG << "number of tracks: " << ntrk << endreq;
00127
00128 for(i=0; i<ntrk; i++){
00129 rectrk = event -> getRecTrk(i);
00130 nhit = rectrk -> getNHits();
00131
00132 log << MSG::DEBUG << "number of hits: " << nhit << endreq;
00133 for(k=0; k<nhit; k++){
00134 rechit = rectrk -> getRecHit(k);
00135 lay = rechit -> getLayid();
00136 cel = rechit -> getCellid();
00137 wir = m_mdcGeomSvc -> Wire(lay, cel) -> Id();
00138 lr = rechit -> getLR();
00139 doca = rechit -> getDocaInc();
00140 dmeas = rechit -> getDmeas();
00141 residual = rechit -> getResiInc();
00142 zhit = rechit -> getZhit();
00143
00144 if((wir<MdcCalTotCell) && (fabs(zhit) < m_zeast[lay])){
00145 bin = (int)(zhit / m_zwid[lay]);
00146 if( 0 == lr ){
00147 m_hl[wir][bin] -> Fill(residual);
00148 }else if( 1 == lr ){
00149 m_hr[wir][bin] -> Fill(residual);
00150 }
00151 }else{
00152 std::cout << "wir: " << wir << std::endl;
00153 }
00154 }
00155 }
00156 return 1;
00157 }
00158
00159 int Wr2dMdcCalib::updateConst(MdcCalibConst* calconst){
00160 IMessageSvc *msgSvc;
00161 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00162 MsgStream log(msgSvc, "Wr2dMdcCalib");
00163 log << MSG::INFO << "Wr2dMdcCalib::updateConst()" << endreq;
00164
00165 int i;
00166 int k;
00167 int bin;
00168 int lay;
00169 int cel;
00170 double dw;
00171
00172
00173 Int_t ierflg;
00174 Int_t istat;
00175 Int_t nvpar;
00176 Int_t nparx;
00177 Double_t fmin;
00178 Double_t edm;
00179 Double_t errdef;
00180 Double_t arglist[10];
00181
00182 TMinuit *gmtrw = new TMinuit(5);
00183 gmtrw -> SetPrintLevel(-1);
00184 gmtrw -> SetFCN(fcnWireParab);
00185 gmtrw -> SetErrorDef(1.0);
00186 gmtrw -> mnparm(0, "xf", 0.0, 0.01, 0, 0, ierflg);
00187 gmtrw -> mnparm(1, "yf", 0.0, 0.01, 0, 0, ierflg);
00188 gmtrw -> mnparm(2, "xb", 0.0, 0.01, 0, 0, ierflg);
00189 gmtrw -> mnparm(3, "yb", 0.0, 0.01, 0, 0, ierflg);
00190 gmtrw -> mnparm(4, "ten", 0.0, 0.1, 0, 0, ierflg);
00191 arglist[0] = 0;
00192 gmtrw -> mnexcm("Set NOW", arglist, 0, ierflg);
00193
00194 double a;
00195 double dx;
00196 double dy;
00197 double dz;
00198 double length;
00199 double ten[] = {15, 15, 15, 16, 16, 17, 17, 18, 14, 14,
00200 19, 19, 24, 24, 31, 31, 37, 37, 45, 45,
00201 46, 47, 47, 47, 47, 48, 48, 48, 48, 49,
00202 49, 49, 49, 50, 50, 50, 51, 51, 51, 52,
00203 52, 52, 52};
00204 double wpar[5];
00205 double wparErr[5];
00206 double sinphiw;
00207 double cosphiw;
00208 double deldw;
00209 double delxw;
00210 double delyw;
00211
00212 int nFit;
00213 Stat_t entryL;
00214 Stat_t entryR;
00215 Double_t hcen;
00216 Double_t cenL[MdcCalWrNBin];
00217 Double_t errL[MdcCalWrNBin];
00218 Double_t cenR[MdcCalWrNBin];
00219 Double_t errR[MdcCalWrNBin];
00220
00221 ofstream fwlog("logWireCor.txt");
00222 ofstream fwpc("wireposCor.txt");
00223 ofstream fwpcErr("wireposErr.txt");
00224
00225 fwpc << setw(5) << "wireId" << setw(15) << "dx_East(mm)"
00226 << setw(15) << "dy_East(mm)" << setw(15) << "dz_East(mm)"
00227 << setw(15) << "dx_West(mm)" << setw(15) << "dy_West(mm)"
00228 << setw(15) << "dz_West(mm)" << endl;
00229
00230 for(i=0; i<MdcCalTotCell; i++){
00231 for(k=0; k<5; k++) wparErr[k] = 999.0;
00232 nFit = 0;
00233 for(bin=0; bin<MdcCalWrNBin; bin++){
00234 entryL = m_hl[i][bin] -> GetEntries();
00235 entryR = m_hr[i][bin] -> GetEntries();
00236 if((entryL < 100) && (entryR < 100)){
00237 fgBIN[bin] = false;
00238 continue;
00239 } else{
00240 fgBIN[bin] = true;
00241 }
00242 nFit++;
00243
00244 if(1 == m_param.fgCalib[lay]){
00245 hcen = m_hl[i][bin] -> GetMean();
00246 if(entryL > 500){
00247 m_hl[i][bin] -> Fit("gaus", "Q", "", hcen-1.5, hcen+1.5);
00248 cenL[bin] = m_hl[i][bin]->GetFunction("gaus")->GetParameter(1);
00249 errL[bin] = m_hl[i][bin]->GetFunction("gaus")->GetParameter(2);
00250 }
00251 else{
00252 cenL[bin] = hcen;
00253 errL[bin] = m_hl[i][bin] -> GetRMS();
00254 }
00255
00256 hcen = m_hr[i][bin] -> GetMean();
00257 if(entryR > 500){
00258 m_hr[i][bin] -> Fit("gaus", "Q", "", hcen-1.5, hcen+1.5);
00259 cenR[bin] = m_hr[i][bin]->GetFunction("gaus")->GetParameter(1);
00260 errR[bin] = m_hr[i][bin]->GetFunction("gaus")->GetParameter(2);
00261 }
00262 else{
00263 cenR[bin] = hcen;
00264 errR[bin] = m_hr[i][bin] -> GetRMS();
00265 }
00266 } else{
00267 cenL[bin] = 0.0;
00268 errL[bin] = 0.15;
00269 cenR[bin] = 0.0;
00270 errR[bin] = 0.15;
00271 }
00272 }
00273 if(nFit < 8) continue;
00274
00275 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
00276 cel = m_mdcGeomSvc -> Wire(i) -> Cell();
00277 zMIN = m_zwest[lay];
00278 zMAX = m_zeast[lay];
00279 wpar[0] = m_mdcGeomSvc->Wire(lay, 0)->Backward().x();
00280 wpar[1] = m_mdcGeomSvc->Wire(lay, 0)->Backward().y();
00281 wpar[2] = m_mdcGeomSvc->Wire(lay, 0)->Forward().x();
00282 wpar[3] = m_mdcGeomSvc->Wire(lay, 0)->Forward().y();
00283 wpar[4] = ten[lay];
00284
00285 a = 9.47e-06 / (2 * wpar[4]);
00286 dx = wpar[0] - wpar[2];
00287 dy = wpar[1] - wpar[3];
00288 dz = zMAX - zMIN;
00289 length = sqrt(dx*dx + dz*dz);
00290
00291 for(bin=0; bin<MdcCalWrNBin; bin++){
00292 zBIN[bin] = m_zbinCen[lay][bin];
00293 xBIN[bin] = (wpar[0]-wpar[2])*(zBIN[bin]-zMIN)/(zMAX-zMIN) + wpar[2];
00294 yBIN[bin] = a*zBIN[bin]*zBIN[bin] + (wpar[1]-wpar[3])*zBIN[bin]/length
00295 + 0.5*(wpar[1]+wpar[3]) - a*length*length/4.0;
00296
00297 sinphiw = yBIN[bin] / sqrt(xBIN[bin]*xBIN[bin] + yBIN[bin]*yBIN[bin]);
00298 cosphiw = xBIN[bin] / sqrt(xBIN[bin]*xBIN[bin] + yBIN[bin]*yBIN[bin]);
00299
00300 deldw = - (cenL[bin]-cenR[bin])/2.0;
00301 delxw = -deldw * sinphiw;
00302 delyw = deldw * cosphiw;
00303
00304 fwlog << setw(3) << lay << setw(4) << cel << setw(5) << i
00305 << setw(4) << bin << setw(15) << zBIN[bin]
00306 << setw(15) << deldw << setw(15) << delxw
00307 << setw(15) << delyw << endl;
00308
00309 xBIN[bin] += delxw;
00310 yBIN[bin] += delyw;
00311
00312 zBINERR[bin] = sqrt((errL[bin]*errL[bin]) + (errR[bin]*errR[bin]))/2.0;
00313 }
00314
00315 arglist[0] = 1;
00316 arglist[1] = wpar[0];
00317 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
00318 arglist[0] = 2;
00319 arglist[1] = wpar[1];
00320 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
00321 arglist[0] = 3;
00322 arglist[1] = wpar[2];
00323 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
00324 arglist[0] = 4;
00325 arglist[1] = wpar[3];
00326 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
00327 arglist[0] = 5;
00328 arglist[1] = wpar[4];
00329 gmtrw -> mnexcm("SET PARameter", arglist, 2, ierflg);
00330 gmtrw -> mnexcm("FIX", arglist, 1, ierflg);
00331
00332 arglist[0] = 2000;
00333 arglist[1] = 0.1;
00334 gmtrw -> mnexcm("MIGRAD", arglist, 2, ierflg);
00335 gmtrw -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
00336
00337 if( (0==ierflg) && (3==istat) ){
00338 gmtrw -> GetParameter(0, wpar[0], wparErr[0]);
00339 gmtrw -> GetParameter(1, wpar[1], wparErr[1]);
00340 gmtrw -> GetParameter(2, wpar[2], wparErr[2]);
00341 gmtrw -> GetParameter(3, wpar[3], wparErr[3]);
00342 gmtrw -> GetParameter(4, wpar[4], wparErr[4]);
00343 }
00344 gmtrw -> mnexcm("RELease", arglist, 1, ierflg);
00345
00346 fwlog << setw(5) << i << setw(15) << wpar[0] << setw(15) << wpar[1]
00347 << setw(15) << wpar[2] << setw(15) << wpar[3] << endl;
00348 fwpc << setw(5) << i << setw(15) << wpar[0] << setw(15) << wpar[1]
00349 << setw(15) << "0" << setw(15) << wpar[2]
00350 << setw(15) << wpar[3] << setw(15) << "0" << endl;
00351 fwpcErr << setw(5) << i << setw(15) << wparErr[0] << setw(15) << wparErr[1]
00352 << setw(15) << wparErr[2] << setw(15) << wparErr[3] << endl;
00353 }
00354 fwlog.close();
00355 fwpc.close();
00356 fwpcErr.close();
00357
00358 delete gmtrw;
00359 return 1;
00360 }
00361
00362 void Wr2dMdcCalib::fcnWireParab(Int_t &npar, Double_t *gin,
00363 Double_t &f, Double_t *par, Int_t iflag){
00364 Int_t bin;
00365 Double_t xfit;
00366 Double_t yfit;
00367
00368 double a = 9.47e-06 / (2 * par[4]);
00369 double dx = par[0] - par[2];
00370 double dy = par[1] - par[3];
00371 double dz = zMAX - zMIN;
00372 double length = sqrt(dx*dx + dz*dz);
00373
00374 Double_t chisq = 0.0;
00375 Double_t deta;
00376
00377 for(bin=0; bin<MdcCalWrNBin; bin++){
00378 if( fgBIN[bin] ){
00379 xfit = (par[0] - par[2]) * (zBIN[bin] - zMIN) / (zMAX - zMIN) + par[2];
00380 yfit = a*zBIN[bin]*zBIN[bin] + (par[1] - par[3])*zBIN[bin]/length
00381 + 0.5*(par[1] + par[3]) - a*length*length/4.0;
00382
00383 deta = ( (xfit-xBIN[bin])*(xfit-xBIN[bin])+
00384 (yfit-yBIN[bin])*(yfit-yBIN[bin]) ) / (zBINERR[bin]*zBINERR[bin]);
00385 chisq += deta;
00386 }
00387 }
00388 f = chisq;
00389 }
00390