BOSS 7.1.0
BESIII Offline Software System
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ResiAlign.cxx
Go to the documentation of this file.
1#include "MdcAlignAlg/ResiAlign.h"
2#include "MdcAlignAlg/MdcAlignAlg.h"
3#include "MdcAlignAlg/Alignment.h"
4
5#include "GaudiKernel/INTupleSvc.h"
6#include "GaudiKernel/IService.h"
7#include "GaudiKernel/IMessageSvc.h"
8#include "GaudiKernel/StatusCode.h"
9#include "GaudiKernel/ISvcLocator.h"
10#include "GaudiKernel/Bootstrap.h"
11
12#include "EventModel/Event.h"
13#include "EventModel/EventHeader.h"
14#include "MdcRawEvent/MdcDigi.h"
15#include "Identifier/Identifier.h"
16#include "Identifier/MdcID.h"
17
18#include "TF1.h"
19#include "TCanvas.h"
20
21#include <iostream>
22#include <fstream>
23#include <iomanip>
24#include <string>
25#include <cstring>
26
27using namespace std;
28
29ResiAlign::ResiAlign(){
30 //m_ndiv = 6;
31 m_ndiv = 12;
32 m_resiCut = 1.0;
33 for(int i=0; i<NEP; i++) m_npoint[i] = 0;
34
35 for(int lay=0; lay<LAYERNMAX; lay++){
36 if(lay < 8){
37 m_docaMin[lay] = 0.6; // mm
38 m_docaMax[lay] = 4.8; // mm
39 } else{
40 m_docaMin[lay] = 0.8; // mm
41 m_docaMax[lay] = 6.4; // mm
42 }
43 }
44 for(int lay=0; lay<LAYERNMAX; lay++){
45 if((0==lay) || (7==lay) || (8==lay) || (19==lay) || (20==lay) ||
46 (35==lay) || (36==lay) || (42==lay) ) m_layBound[lay] = true;
47 else m_layBound[lay] = false;
48 }
49
50 m_ncut1 = 0;
51 m_ncut2 = 0;
52 m_ncut3 = 0;
53 m_ncut4 = 0;
54 m_ncut5 = 0;
55 m_ncut6 = 0;
56 m_ncut7 = 0;
57 m_ncut8 = 0;
58 m_ncut9 = 0;
59 m_ncut10 = 0;
60 m_ncut11 = 0;
61 m_ncut12 = 0;
62 m_ncut13 = 0;
63}
64
65ResiAlign::~ResiAlign(){
66}
67
68void ResiAlign::clear(){
69 delete m_hresAll;
70 delete m_hresInn;
71 delete m_hresStp;
72 delete m_hresOut;
73 for(int lay=0; lay<LAYERNMAX; lay++){
74 delete m_hresLay[lay];
75 for(int i=0; i<4; i++) delete m_hresLay_LR[lay][i];
76 }
77 for(int i=0; i<NEP; i++) delete m_gr[i];
78}
79
80void ResiAlign::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
81 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc){
82 IMessageSvc* msgSvc;
83 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
84 MsgStream log(msgSvc, "ResiAlign");
85 log << MSG::INFO << "ResiAlign::initialize()" << endreq;
86
87 m_hlist = hlist;
88 m_mdcGeomSvc = mdcGeomSvc;
89 m_mdcFunSvc = mdcFunSvc;
90 m_mdcUtilitySvc = mdcUtilitySvc;
91
92 double zeast;
93 for(int lay=0; lay<43; lay++){
94 zeast = m_mdcGeomSvc->Wire(lay, 0)->Backward().z();
95 m_zrange[lay][1] = 2.0 * fabs(zeast) / (double)m_ndiv;
96 m_zrange[lay][0] = -1.0 * m_zrange[lay][1];
97
98 m_radii[lay] = m_mdcGeomSvc->Layer(lay)->Radius();
99 }
100
101 for(int wir=0; wir<WIRENMAX; wir++){
102 m_xe[wir] = m_mdcGeomSvc->Wire(wir)->Backward().x();
103 m_ye[wir] = m_mdcGeomSvc->Wire(wir)->Backward().y();
104 m_ze[wir] = m_mdcGeomSvc->Wire(wir)->Backward().z();
105 m_xw[wir] = m_mdcGeomSvc->Wire(wir)->Forward().x();
106 m_yw[wir] = m_mdcGeomSvc->Wire(wir)->Forward().y();
107 m_zw[wir] = m_mdcGeomSvc->Wire(wir)->Forward().z();
108 }
109
110 char hname[200];
111 int iEP;
112
113 INTupleSvc* ntupleSvc;
114 Gaudi::svcLocator() -> service("NTupleSvc", ntupleSvc);
115 for(iEP=0; iEP<=NEP; iEP++){
116 if(iEP < NEP) sprintf(hname, "FILE137/align%02d", iEP);
117 else sprintf(hname, "FILE137/alignAll");
118
119 NTuplePtr nt(ntupleSvc, hname);
120 if( nt ) m_tuple[iEP] = nt;
121 else{
122 m_tuple[iEP] = ntupleSvc->book(hname, CLID_ColumnWiseTuple,"align");
123 if (m_tuple[iEP]) {
124 m_tuple[iEP]->addItem ("run", m_iRun[iEP]);
125 m_tuple[iEP]->addItem ("evt", m_iEvt[iEP]);
126 m_tuple[iEP]->addItem ("resi", m_resi[iEP]);
127 m_tuple[iEP]->addItem ("p", m_p[iEP]);
128 m_tuple[iEP]->addItem ("pt", m_pt[iEP]);
129 m_tuple[iEP]->addItem ("phi", m_phi[iEP]);
130 m_tuple[iEP]->addItem ("lay", m_lay[iEP]);
131 m_tuple[iEP]->addItem ("lr", m_lr[iEP]);
132 m_tuple[iEP]->addItem ("cel", m_cel[iEP]);
133 }
134 else {
135 log << MSG::FATAL << "Cannot book N-tuple:"
136 << long(m_tuple[iEP]) << endmsg;
137 }
138 }
139 }
140
141 m_hnTrk = new TH1F("HNtrack", "", 10, -0.5, 9.5);
142 m_hlist->Add(m_hnTrk);
143
144 m_hnHit = new TH1F("HNhit", "", 100, -0.5, 99.5);
145 m_hlist->Add(m_hnHit);
146
147 m_hlayHitmap = new TH1F("Hitmap", "", 43, -0.5, 42.5);
148 m_hlist->Add(m_hlayHitmap);
149
150 m_hresAll = new TH1F("HResAllInc", "", 200, -1.0, 1.0);
151 m_hlist->Add(m_hresAll);
152
153 m_hresInn = new TH1F("HResInnInc", "", 200, -1.0, 1.0);
154 m_hlist->Add(m_hresInn);
155
156 m_hresStp = new TH1F("HResStpInc", "", 200, -1.0, 1.0);
157 m_hlist->Add(m_hresStp);
158
159 m_hresOut = new TH1F("HResOutInc", "", 200, -1.0, 1.0);
160 m_hlist->Add(m_hresOut);
161
162 int lay;
163 for(lay=0; lay<LAYERNMAX; lay++){
164 sprintf(hname, "Res_Layer%02d", lay);
165 m_hresLay[lay] = new TH1F(hname, "", 200, -1.0, 1.0);
166 m_hlist->Add(m_hresLay[lay]);
167
168 for(int i=0; i<4; i++){
169 if(0==i) sprintf(hname, "Resi_Lay%02d_Up_L", lay);
170 else if(1==i) sprintf(hname, "Resi_Lay%02d_Up_R", lay);
171 else if(2==i) sprintf(hname, "Resi_Lay%02d_Dw_L", lay);
172 else sprintf(hname, "Resi_Lay%02d_Dw_R", lay);
173 m_hresLay_LR[lay][i] = new TH1F(hname, "", 200, -1.0, 1.0);
174 m_hlist->Add(m_hresLay_LR[lay][i]);
175 }
176 }
177
178 for(iEP=0; iEP<NEP; iEP++){
179 m_gr[iEP] = new TGraph();
180 sprintf(hname, "grResi%02d", iEP);
181 m_gr[iEP]->SetName(hname);
182 m_hlist->Add(m_gr[iEP]);
183 }
184 m_fevt.open("evt.txt");
185}
186
187bool ResiAlign::fillHist(MdcAliEvent* event){
188 IMessageSvc* msgSvc;
189 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
190 MsgStream log(msgSvc, "ResiAlign");
191 log << MSG::DEBUG << "ResiAlign::fillHist()" << endreq;
192
193 bool esCutFg = event->getEsCutFlag();
194 if( ! esCutFg ){
195 m_ncut1++;
196 return true;
197 }
198
199 int i = 0;
200 int k;
201
202 int trkStat;
203 double dr;
204 double phi0;
205 double kappa;
206 double dz;
207 double tanl;
208 double chisq;
209 double p;
210 double pt;
211
212 int nhits;
213 int lay;
214 int cel;
215 int wir;
216 int lr;
217 int iEnd;
218 int iEP;
219
220 double doca;
221 double resi;
222 double zhit;
223 double wphi;
224 double dphi;
225 double hitPhi;
226 double xx;
227 double yy;
228 double rr;
229 int stat;
230 MdcAliRecTrk* rectrk;
231 MdcAliRecHit* rechit;
232 int nhitlay;
233 bool fgHitLay[LAYERNMAX];
234
235 IDataProviderSvc* eventSvc = NULL;
236 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
237 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
238 if (!eventHeader) {
239 log << MSG::FATAL << "Could not find Event Header" << endreq;
240 m_ncut3++;
241 return( StatusCode::FAILURE);
242 }
243 int iEvt = eventHeader->eventNumber();
244 int iRun = eventHeader->runNumber();
245
246 int nTrk = event -> getNTrk();
247 m_hnTrk->Fill(nTrk);
248 m_fevt << setw(10) << iRun << setw(10) << iEvt << setw(10) << nTrk << endl;
249 if((nTrk < m_param.nTrkCut[0]) || (nTrk > m_param.nTrkCut[1])){
250 m_ncut2++;
251 return true;
252 }
253
254 for(i=0; i<nTrk; i++){
255 rectrk = event->getRecTrk(i);
256 nhits = rectrk->getNHits();
257 trkStat = rectrk->getStat();
258// if (0 != trkStat) continue;
259
260 // dr cut
261 dr = rectrk -> getDr();
262 if(fabs(dr) > m_param.drCut){ m_ncut4++; continue; }
263
264 phi0 = rectrk -> getPhi0();
265 kappa = rectrk -> getKappa();
266
267 // dz cut
268 dz = rectrk -> getDz();
269 if(fabs(dz) > m_param.dzCut){ m_ncut5++; continue; }
270
271 for(lay=0; lay<LAYERNMAX; lay++){
272 fgHitLay[lay] = false;
273 }
274
275 m_hnHit->Fill(nhits);
276 for(k=0; k<nhits; k++){
277 rechit = rectrk -> getRecHit(k);
278 lay = rechit -> getLayid();
279 fgHitLay[lay] = true;
280 }
281
282 nhitlay = 0;
283 for(lay=0; lay<LAYERNMAX; lay++){
284 if(fgHitLay[lay]) nhitlay++;
285 }
286// cout << "hitlay " << nhitlay << endl;
287 if(nhitlay < m_param.nHitLayCut){ m_ncut6++; continue; }
288
289 tanl = rectrk -> getTanLamda();
290 chisq = rectrk -> getChisq();
291 p = rectrk -> getP();
292 pt = rectrk -> getPt();
293
294 if((fabs(pt)<m_param.ptCut[0]) || (fabs(pt)>m_param.ptCut[1])){ m_ncut7++; continue;}
295
296 HepVector helix = rectrk->getHelix();
297 for(k=0; k<nhits; k++){
298 rechit = rectrk->getRecHit(k);
299 lay = rechit->getLayid();
300 cel = rechit->getCellid();
301 lr = rechit->getLR();
302 doca = rechit -> getDocaInc();
303 zhit = rechit->getZhit();
304
305 stat = rechit -> getStat();
306 if((1 == m_param.hitStatCut) && (1 != stat)){ m_ncut8++; continue; }
307
308 if (1 == m_param.resiType) {
309 resi = rechit->getResiExcLR();
310 } else {
311 resi = rechit->getResiIncLR();
312 }
313 resi *= -1.0;
314 //maqm if( (1==isnan(resi)) || (fabs(resi) > m_resiCut) ||
315 if( (1==std::isnan(resi)) || (fabs(resi) > m_resiCut) ||
316 (fabs(doca) > m_docaMax[lay]) || (fabs(doca) < m_docaMin[lay]) ){
317 m_ncut9++;
318 continue;
319 }
320
321 if(m_param.fgAdjacLayerCut){
322 if(0 == lay){
323 if( ! fgHitLay[1] ){ m_ncut10++; continue; }
324 } else if(42 == lay){
325 if( ! fgHitLay[41] ){ m_ncut11++; continue; }
326 } else{
327 if( (!fgHitLay[lay-1]) && (!fgHitLay[lay+1]) ){ m_ncut12++; continue; }
328 // for boundary layers
329 if( m_param.fgBoundLayerCut && m_layBound[lay] &&
330 ((!fgHitLay[lay-1]) || (!fgHitLay[lay+1])) ){
331 m_ncut13++;
332 continue;
333 }
334 }
335 }
336 m_hlayHitmap->Fill(lay);
337
338 // fill alignment trees
339 if((zhit < m_zrange[lay][0]) || (zhit > m_zrange[lay][1])){
340 wir = m_mdcGeomSvc -> Wire(lay, cel) -> Id();
341 xx = (zhit - m_zw[wir]) * (m_xe[wir] - m_xw[wir]) /
342 (m_ze[wir] - m_zw[wir]) + m_xw[wir];
343 yy = (zhit - m_zw[wir]) * (m_ye[wir] - m_yw[wir]) /
344 (m_ze[wir] - m_zw[wir]) + m_yw[wir];
345 rr = sqrt( (xx * xx) + (yy * yy) );
346 dphi = fabs(doca) / m_radii[lay];
347
348 if( yy >= 0 ) wphi = acos(xx / rr);
349 else wphi = PI2 - acos(xx / rr);
350 if(1 == lr) hitPhi = wphi + dphi;
351 else hitPhi = wphi - dphi;
352 if(hitPhi < 0) hitPhi += PI2;
353 else if(hitPhi > PI2) hitPhi -= PI2;
354
355 if(zhit < m_zrange[lay][0]) iEnd = 0; // west
356 else iEnd = 1; // east
357 iEP = Alignment::getEpId(lay, iEnd);
358
359 m_iRun[iEP] = iRun;
360 m_iEvt[iEP] = iEvt;
361 m_resi[iEP] = resi;
362 m_p[iEP] = p;
363 m_pt[iEP] = pt;
364 m_phi[iEP] = hitPhi;
365 m_lay[iEP] = lay;
366 m_lr[iEP] = lr;
367 m_cel[iEP] = cel;
368 m_tuple[iEP]->write();
369
370 m_resi[NEP] = resi;
371 m_p[NEP] = p;
372 m_pt[NEP] = pt;
373 m_phi[NEP] = hitPhi;
374 m_lay[NEP] = lay;
375 m_lr[NEP] = lr;
376 m_cel[NEP] = cel;
377 m_tuple[NEP]->write();
378
379 m_hresAll->Fill(resi);
380 if(lay < 8) m_hresInn->Fill(resi);
381 else if(lay < 20) m_hresStp->Fill(resi);
382 else m_hresOut->Fill(resi);
383 m_hresLay[lay]->Fill(resi);
384 if(0==lr){
385 if(hitPhi<3.14) m_hresLay_LR[lay][0]->Fill(resi); // up L
386 else m_hresLay_LR[lay][2]->Fill(resi); // down L
387 } else{
388 if(hitPhi<3.14) m_hresLay_LR[lay][1]->Fill(resi); // up R
389 else m_hresLay_LR[lay][3]->Fill(resi); // down R
390 }
391
392 m_gr[iEP]->SetPoint(m_npoint[iEP], hitPhi, resi);
393 m_npoint[iEP]++;
394 }
395 }
396 }
397
398 return true;
399}
400
401void ResiAlign::updateConst(MdcAlignPar* alignPar){
402 IMessageSvc* msgSvc;
403 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
404 MsgStream log(msgSvc, "ResiAlign");
405 log << MSG::INFO << "ResiAlign::updateConst()" << endreq;
406 m_fevt.close();
407
408 int iEP;
409 double par[3];
410 double err[3];
411 double dx;
412 double dy;
413 double rz;
414 double rLayer[] = { 120.225, 205.0, 237.55, 270.175,
415 302.625, 334.775, 366.65, 500.0,
416 120.225, 205.0, 237.55, 270.175,
417 302.625, 334.775, 366.65, 500.0 };
418
419 TCanvas c1("c1", "c1", 10, 10, 700, 500);
420
421 TF1* fResPhi = new TF1("fResPhi", funResi, 0, PI2, 3);
422 fResPhi->SetParameter(0, 0.0);
423 fResPhi->SetParameter(1, 0.0);
424 fResPhi->SetParameter(2, 0.0);
425
426 for(iEP=0; iEP<NEP; iEP++){
427 if((m_gr[iEP]->GetN()) > 500){
428 m_gr[iEP]->Fit("fResPhi", "V");
429 par[0] = fResPhi->GetParameter(0);
430 par[1] = fResPhi->GetParameter(1);
431 par[2] = fResPhi->GetParameter(2);
432
433 err[0] = fResPhi->GetParError(0);
434 err[1] = fResPhi->GetParError(1);
435 err[2] = fResPhi->GetParError(2);
436
437 dx = -1.0 * par[1];
438 dy = par[2];
439 rz = par[0] / rLayer[iEP];
440
441 // assume the shift of the outer section is 0
442 if (7==iEP || 15==iEP) {
443 dx = 0.0;
444 dy = 0.0;
445 rz = 0.0;
446 par[0] = 0.0;
447 par[1] = 0.0;
448 par[2] = 0.0;
449 }
450 alignPar->setDelDx(iEP, dx);
451 alignPar->setDelDy(iEP, dy);
452 alignPar->setDelRz(iEP, rz);
453
454 alignPar->setErrDx(iEP, err[1]);
455 alignPar->setErrDy(iEP, err[2]);
456 alignPar->setErrRz(iEP, err[0]/rLayer[iEP]);
457 }
458 }
459
460 cout << "TrackCut: cut1: " << m_ncut1 << ", cut2: " << m_ncut2 << ", cut3: " << m_ncut3
461 << ", cut4: " << m_ncut4 << ", cut5: " << m_ncut5 << ", cut6: " << m_ncut6
462 << ", cut7: " << m_ncut7 << endl;
463 cout << "HitCut: cut8: " << m_ncut8 << ", cut9: " << m_ncut9 << ", cut10: " << m_ncut10
464 << ", cut11: " << m_ncut11 << ", cut12: " << m_ncut12 << ", cut13: " << m_ncut13 << endl;
465
466 delete fResPhi;
467}
468
469Double_t ResiAlign::funResi(Double_t* x, Double_t* par){
470 Double_t val;
471 val = par[0] + par[1]*sin(x[0]) + par[2]*cos(x[0]);
472 return val;
473}
474
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
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)
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
NULL
#define NULL
Definition: VertexExtrapolate.h:14
IMdcGeomSvc::Layer
virtual const MdcGeoLayer *const Layer(unsigned id)=0
IMdcGeomSvc::Wire
virtual const MdcGeoWire *const Wire(unsigned id)=0
MdcAliParams::dzCut
double dzCut
Definition: MdcAliParams.h:29
MdcAliParams::drCut
double drCut
Definition: MdcAliParams.h:28
MdcAliParams::nTrkCut
int nTrkCut[2]
Definition: MdcAliParams.h:21
MdcAliParams::fgAdjacLayerCut
int fgAdjacLayerCut
Definition: MdcAliParams.h:19
MdcAliParams::ptCut
double ptCut[2]
Definition: MdcAliParams.h:26
MdcAliParams::fgBoundLayerCut
int fgBoundLayerCut
Definition: MdcAliParams.h:20
MdcAliRecHit::getCellid
int getCellid() const
Definition: MdcAliRecHit.h:21
MdcAliRecHit::getLR
int getLR() const
Definition: MdcAliRecHit.h:22
MdcAliRecHit::getResiExcLR
double getResiExcLR() const
Definition: MdcAliRecHit.h:34
MdcAliRecHit::getLayid
int getLayid() const
Definition: MdcAliRecHit.h:20
MdcAliRecHit::getZhit
double getZhit() const
Definition: MdcAliRecHit.h:39
MdcAliRecHit::getResiIncLR
double getResiIncLR() const
Definition: MdcAliRecHit.h:31
MdcAliRecTrk::getRecHit
MdcAliRecHit * getRecHit(int index) const
Definition: MdcAliRecTrk.h:38
MdcAliRecTrk::getStat
int getStat() const
Definition: MdcAliRecTrk.h:26
MdcAliRecTrk::getNHits
int getNHits() const
Definition: MdcAliRecTrk.h:37
MdcAliRecTrk::getHelix
HepVector getHelix() const
Definition: MdcAliRecTrk.h:32
MdcAlignPar::setErrRz
void setErrRz(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:40
MdcAlignPar::setDelDy
void setDelDy(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:29
MdcAlignPar::setDelRz
void setDelRz(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:33
MdcAlignPar::setErrDx
void setErrDx(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:35
MdcAlignPar::setDelDx
void setDelDx(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:28
MdcAlignPar::setErrDy
void setErrDy(int iEP, double val)
Definition: MdcAlignAlg/MdcAlignPar.h:36
MdcGeoLayer::Radius
double Radius(void) const
Definition: MdcGeoLayer.h:160
MdcGeoWire::Forward
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
MdcGeoWire::Backward
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128
ResiAlign::funResi
static Double_t funResi(double *x, double *par)
ResiAlign::initialize
void initialize(TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)
Definition: ResiAlign.cxx:80
ResiAlign::fillHist
bool fillHist(MdcAliEvent *event)
Definition: ResiAlign.cxx:187
ResiAlign::clear
void clear()
Definition: ResiAlign.cxx:68
ResiAlign::updateConst
void updateConst(MdcAlignPar *alignPar)
Definition: ResiAlign.cxx:401
Alignment::PI2
const double PI2
Definition: Alignment.h:41
Alignment::LAYERNMAX
const int LAYERNMAX
Definition: Alignment.h:45
Alignment::NEP
const int NEP
Definition: Alignment.h:48
Alignment::getEpId
int getEpId(int lay, int iEnd)
Definition: Alignment.cxx:18
Alignment::WIRENMAX
const int WIRENMAX
Definition: Alignment.h:44
nhits
int nhits
Definition: prepare_data.cxx:29