BOSS 7.1.0
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
159void Wr2dMdcCalib::printCut() const{
160 MdcCalib::printCut();
161}
162
163int Wr2dMdcCalib::updateConst(MdcCalibConst* calconst){
164 IMessageSvc *msgSvc;
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;
171 int bin;
172 int lay;
173 int cel;
174 double dw;
175
176 // minuit for wires
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);
188 gmtrw -> SetFCN(fcnWireParab);
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;
220 Double_t cenL[MdcCalWrNBin];
221 Double_t errL[MdcCalWrNBin];
222 Double_t cenR[MdcCalWrNBin];
223 Double_t errR[MdcCalWrNBin];
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
234 for(i=0; i<MdcCalTotCell; i++){
235 for(k=0; k<5; k++) wparErr[k] = 999.0;
236 nFit = 0;
237 for(bin=0; bin<MdcCalWrNBin; bin++){
238 entryL = m_hl[i][bin] -> GetEntries();
239 entryR = m_hr[i][bin] -> GetEntries();
240 if((entryL < 100) && (entryR < 100)){
241 fgBIN[bin] = false;
242 continue;
243 } else{
244 fgBIN[bin] = true;
245 }
246 nFit++;
247
248 if(1 == m_param.fgCalib[lay]){
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{
256 cenL[bin] = hcen;
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{
267 cenR[bin] = hcen;
268 errR[bin] = m_hr[i][bin] -> GetRMS();
269 }
270 } else{
271 cenL[bin] = 0.0;
272 errL[bin] = 0.15;
273 cenR[bin] = 0.0;
274 errR[bin] = 0.15;
275 }
276 }
277 if(nFit < 8) continue;
278
279 lay = m_mdcGeomSvc -> Wire(i) -> Layer();
280 cel = m_mdcGeomSvc -> Wire(i) -> Cell();
281 zMIN = m_zwest[lay];
282 zMAX = m_zeast[lay];
283 wpar[0] = m_mdcGeomSvc->Wire(lay, 0)->Backward().x();
284 wpar[1] = m_mdcGeomSvc->Wire(lay, 0)->Backward().y();
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];
292 dz = zMAX - zMIN;
293 length = sqrt(dx*dx + dz*dz);
294
295 for(bin=0; bin<MdcCalWrNBin; bin++){
296 zBIN[bin] = m_zbinCen[lay][bin];
297 xBIN[bin] = (wpar[0]-wpar[2])*(zBIN[bin]-zMIN)/(zMAX-zMIN) + wpar[2];
298 yBIN[bin] = a*zBIN[bin]*zBIN[bin] + (wpar[1]-wpar[3])*zBIN[bin]/length
299 + 0.5*(wpar[1]+wpar[3]) - a*length*length/4.0;
300
301 sinphiw = yBIN[bin] / sqrt(xBIN[bin]*xBIN[bin] + yBIN[bin]*yBIN[bin]);
302 cosphiw = xBIN[bin] / sqrt(xBIN[bin]*xBIN[bin] + yBIN[bin]*yBIN[bin]);
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
309 << setw(4) << bin << setw(15) << zBIN[bin]
310 << setw(15) << deldw << setw(15) << delxw
311 << setw(15) << delyw << endl;
312
313 xBIN[bin] += delxw;
314 yBIN[bin] += delyw;
315
316 zBINERR[bin] = sqrt((errL[bin]*errL[bin]) + (errR[bin]*errR[bin]))/2.0;
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}
365
366void Wr2dMdcCalib::fcnWireParab(Int_t &npar, Double_t *gin,
367 Double_t &f, Double_t *par, Int_t iflag){
368 Int_t bin;
369 Double_t xfit;
370 Double_t yfit;
371
372 double a = 9.47e-06 / (2 * par[4]);
373 double dx = par[0] - par[2];
374 double dy = par[1] - par[3];
375 double dz = zMAX - zMIN;
376 double length = sqrt(dx*dx + dz*dz);
377
378 Double_t chisq = 0.0;
379 Double_t deta;
380
381 for(bin=0; bin<MdcCalWrNBin; bin++){
382 if( fgBIN[bin] ){
383 xfit = (par[0] - par[2]) * (zBIN[bin] - zMIN) / (zMAX - zMIN) + par[2];
384 yfit = a*zBIN[bin]*zBIN[bin] + (par[1] - par[3])*zBIN[bin]/length
385 + 0.5*(par[1] + par[3]) - a*length*length/4.0;
386
387 deta = ( (xfit-xBIN[bin])*(xfit-xBIN[bin])+
388 (yfit-yBIN[bin])*(yfit-yBIN[bin]) ) / (zBINERR[bin]*zBINERR[bin]);
389 chisq += deta;
390 }
391 }
392 f = chisq;
393}
394
Fill
curve Fill()
sprintf
sprintf(cut,"kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
Fit
void Fit()
Definition: Eangle1D/Fit.cxx:3
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
MdcCalNLayer
const int MdcCalNLayer
Definition: MdcCalParams.h:6
MdcCalTotCell
const int MdcCalTotCell
Definition: MdcCalParams.h:9
MdcCalWrNBin
const int MdcCalWrNBin
Definition: MdcCalParams.h:22
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
IMdcGeomSvc::Wire
virtual const MdcGeoWire *const Wire(unsigned id)=0
MdcCalParams::fgCalib
int fgCalib[MdcCalNLayer]
Definition: MdcCalParams.h:80
MdcCalib::printCut
virtual void printCut() const =0
Definition: MdcCalib.cxx:1407
MdcCalib::clear
virtual void clear()=0
Definition: MdcCalib.cxx:80
MdcCalib::initialize
virtual void initialize(TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)=0
Definition: MdcCalib.cxx:214
MdcCalib::fillHist
virtual int fillHist(MdcCalEvent *event)=0
Definition: MdcCalib.cxx:730
MdcGeoWire::Forward
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
MdcGeoWire::Backward
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128
Wr2dMdcCalib::zBIN
static double zBIN[MdcCalWrNBin]
Definition: Wr2dMdcCalib.h:25
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:366
Wr2dMdcCalib::initialize
void initialize(TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)
Definition: Wr2dMdcCalib.cxx:46
Wr2dMdcCalib::printCut
void printCut() const
Definition: Wr2dMdcCalib.cxx:159
Wr2dMdcCalib::xBIN
static double xBIN[MdcCalWrNBin]
Definition: Wr2dMdcCalib.h:23
Wr2dMdcCalib::zBINERR
static double zBINERR[MdcCalWrNBin]
Definition: Wr2dMdcCalib.h:26
Wr2dMdcCalib::fgBIN
static bool fgBIN[MdcCalWrNBin]
Definition: Wr2dMdcCalib.h:22
Wr2dMdcCalib::yBIN
static double yBIN[MdcCalWrNBin]
Definition: Wr2dMdcCalib.h:24
Wr2dMdcCalib::updateConst
int updateConst(MdcCalibConst *calconst)
Definition: Wr2dMdcCalib.cxx:163
Wr2dMdcCalib::zMAX
static double zMAX
Definition: Wr2dMdcCalib.h:28
Wr2dMdcCalib::zMIN
static double zMIN
Definition: Wr2dMdcCalib.h:27