CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcGeomSvc.cxx
Go to the documentation of this file.
1#include "GaudiKernel/Algorithm.h"
3#include "GaudiKernel/Kernel.h"
4#include "GaudiKernel/IIncidentSvc.h"
5#include "GaudiKernel/Incident.h"
6#include "GaudiKernel/IIncidentListener.h"
7#include "GaudiKernel/IInterface.h"
8#include "GaudiKernel/StatusCode.h"
9#include "GaudiKernel/SvcFactory.h"
10#include "GaudiKernel/PropertyMgr.h"
11#include "GaudiKernel/SmartDataPtr.h"
14#include "GaudiKernel/MsgStream.h"
16#include "EventModel/Event.h"
17#include "GaudiKernel/ISvcLocator.h"
18#include "GaudiKernel/IDataProviderSvc.h"
19#include "GaudiKernel/Bootstrap.h"
20#include <fstream>
21#include <float.h>
22
23
24
25// initialize static members
26bool MdcGeomSvc::m_doSag = true;
29
30MdcGeomSvc::MdcGeomSvc( const std::string& name, ISvcLocator* svcloc ) : Service(name, svcloc) {
31 if(getenv("MDCGEOMSVCROOT")){
32 m_alignFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/MdcAlignPar.dat");
33 //std::cout<<" the MDC alignment file: "<<m_alignFilePath<<std::endl;
34
35 m_wirePosFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/WirePosCalib.dat");
36 //std::cout<<" the MDC wire position file: "<<m_wirePosFilePath<<std::endl;
37
38 m_wireTensionFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/mdcWireTension.dat");
39 //std::cout<<" the MDC wire tension file: "<<m_wireTensionFilePath<<std::endl;
40
41 }
42 else {
43 std::cout<<"A fatal error, contact wangjk..."<<std::endl;
44 }
45
46 declareProperty("doSag", m_doSag = true);
47 declareProperty("readAlignParDataBase", m_readAlignParDataBase = true);
48 declareProperty("mcnoalignment",m_nomcalignment=true);
49 declareProperty("wholeShiftX",m_wholeShiftX = 0.);
50 declareProperty("wholeShiftY",m_wholeShiftY = 0.);
51 declareProperty("wholeShiftZ",m_wholeShiftZ = 0.);
52 declareProperty("wholeRotatX",m_wholeRotatX = 0.);
53 declareProperty("wholeRotatY",m_wholeRotatY = 0.);
54 declareProperty("wholeRotatZ",m_wholeRotatZ = 0.);
55 declareProperty("alignFilePath",m_alignFilePath);
56 declareProperty("wirePosFilePath",m_wirePosFilePath);
57 declareProperty("wireTensionFilePath",m_wireTensionFilePath);
58 declareProperty("useCGEM",m_useCgem = true);
59
60
61}
62
63StatusCode MdcGeomSvc::queryInterface (const InterfaceID& riid, void** ppvInterface ){
64
65 if ( IID_IMdcGeomSvc.versionMatch(riid) ) {
66 *ppvInterface = static_cast<IMdcGeomSvc*> (this);
67 } else {
68 return Service::queryInterface(riid, ppvInterface) ;
69 }
70 return StatusCode::SUCCESS;
71}
72
73StatusCode MdcGeomSvc::initialize ( ) {
74 MsgStream log(messageService(), name());
75 log << MSG::INFO << name() << ": Start of run initialisation" << endreq;
76
77 StatusCode sc = Service::initialize();
78 if ( sc.isFailure() ) return sc;
79 m_mindex=0;
80 m_updataalign = false;
81 IIncidentSvc* incsvc;
82 sc = service("IncidentSvc", incsvc);
83 int priority = 100;
84 if( sc.isSuccess() ){
85 incsvc -> addListener(this, "NewRun", priority);
86 }
87
88 // ReadFilePar(); // get geometry data from file SimUtil/dat/Mdc.txt
89 // Fill(); // get geometry data from Database
90
91 sc = service("CalibDataSvc", m_pCalibDataSvc, true);
92
93 if ( !sc.isSuccess() ) {
94 log << MSG::ERROR
95 << "Could not get IDataProviderSvc interface of CalibXmlCnvSvc"
96 << endreq;
97 return sc;
98 } else {
99 log << MSG::DEBUG
100 << "Retrieved IDataProviderSvc interface of CalibXmlCnvSvc"
101 << endreq;
102 }
103 ReadFilePar();
104 return StatusCode::SUCCESS;
105}
106
107StatusCode MdcGeomSvc::finalize ( ) {
108 MsgStream log(messageService(), name());
109 log << MSG::INFO << name() << ": End of Run" << endreq;
110 return StatusCode::SUCCESS;
111}
112
114 for(vector<MdcGeoLayer*>::iterator it1 = fLayers.begin(); it1 != fLayers.end(); it1++) delete *it1;
115 for(vector<MdcGeoSuper*>::iterator it2 = fSupers.begin(); it2 != fSupers.end(); it2++) delete *it2;
116 for(vector<MdcGeoWire*>::iterator it3 = fWires.begin(); it3 != fWires.end(); it3++) delete *it3;
117 for(vector<MdcGeoEnd*>::iterator it4 = fEnd.begin(); it4 != fEnd.end(); it4++) delete *it4;
118 fGenerals.clear();
119 fWires.clear();
120 fLayers.clear();
121 fSupers.clear();
122 fEnd.clear();
123}
124
125void MdcGeomSvc::clean(){
126 for(vector<MdcGeoLayer*>::iterator it1 = fLayers.begin(); it1 != fLayers.end(); it1++) delete *it1;
127 for(vector<MdcGeoSuper*>::iterator it2 = fSupers.begin(); it2 != fSupers.end(); it2++) delete *it2;
128 for(vector<MdcGeoWire*>::iterator it3 = fWires.begin(); it3 != fWires.end(); it3++) delete *it3;
129 for(vector<MdcGeoEnd*>::iterator it4 = fEnd.begin(); it4 != fEnd.end(); it4++) delete *it4;
130 fGenerals.clear();
131 fWires.clear();
132 fLayers.clear();
133 fSupers.clear();
134 fEnd.clear();
135}
136
137
138void MdcGeomSvc::ReadFilePar(){
139 std::string geometryFilePath = getenv("MDCSIMROOT");
140 if(m_useCgem) geometryFilePath += "/dat/Mdc_cgem.txt";
141 else geometryFilePath += "/dat/Mdc.txt";
142 std::ifstream inFile(geometryFilePath.c_str() );
143
144 if(!inFile.good()){
145 std::cout<<"Error, mdc parameters file not exist"<<std::endl;
146 return;
147 }
148 std::string alignFilePath;
149 std::string wirePosFilePath;
150 std::string wireTensionFilePath;
151 if (!m_updataalign) {
152 alignFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/MdcAlignPar.dat");
153 wirePosFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/WirePosCalib.dat");
154 wireTensionFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/mdcWireTension.dat");
155 } else {
156 alignFilePath = m_alignFilePath;
157 wirePosFilePath = m_wirePosFilePath;
158 wireTensionFilePath = m_wireTensionFilePath;
159 }
160 std::ifstream alignFile(alignFilePath.c_str());
161 std::ifstream wirePosFile(wirePosFilePath.c_str());
162 std::ifstream wireTensionFile(wireTensionFilePath.c_str());
163
164 if(!wirePosFile){
165 std::cout<<"Error, mdc wire position file not exist..."<<std::endl;
166 }
167
168 if(!wireTensionFile){
169 std::cout<<"Error, mdc wire tension file not exist..."<<std::endl;
170 }
171
172 std::vector<double> wireTensionVec;
173 //cout << __FILE__ << " " << __LINE__ << " readdb = " << m_readAlignParDataBase
174 // << " update align = " << m_updataalign << endl;
175 if (m_readAlignParDataBase && m_updataalign) {
176 ReadTensionDataBase(wireTensionVec);
177 }else {
178 cout << "Read wireTension from file:" << wireTensionFilePath.c_str() << endl;
179 int idd;
180 double tension;
181 for(int ii=0; ii<6796; ii++){
182 wireTensionFile>>idd>>tension;
183
184 if(getSagFlag()){
185 wireTensionVec.push_back(tension);
186 }
187 else {
188 tension = DBL_MAX;
189 wireTensionVec.push_back(tension);
190 }
191 }
192 }
193
194 std::vector<vector<double> > wirePosVec(6796);
195 if (m_readAlignParDataBase && m_updataalign) {
196 ReadWirePosDataBase(wirePosVec);
197 } else {
198 cout << "Read wirePosition from file:" << wirePosFilePath.c_str() << endl;
199 double dxe,dye,dze,dxw,dyw,dzw;
200 int wid;
201 std::string title;
202 getline(wirePosFile, title);
203 wirePosFile.seekg(1,ios::cur);
204 for(int j=0; j<6796; j++){
205 wirePosFile>>wid>>dxe>>dye>>dze>>dxw>>dyw>>dzw;
206 wirePosVec[j].push_back(dxe);
207 wirePosVec[j].push_back(dye);
208 wirePosVec[j].push_back(dze);
209 wirePosVec[j].push_back(dxw);
210 wirePosVec[j].push_back(dyw);
211 wirePosVec[j].push_back(dzw);
212 /// test
213 /*
214 for(int jj=0;jj<6;jj++){
215 std::cout<<" wid: "<<wid<<" dxe: "<<wirePosVec[j][0]<<" dye: "<<wirePosVec[j][1]<<" dze: "<<wirePosVec[j][2]<<" dxw: "<<wirePosVec[j][3]<<" dyw: "<<wirePosVec[j][4]<<" dzw: "<<wirePosVec[j][5]<<std::endl;
216 }
217 */
218 }
219 }
220
221 double signalWireR, fieldWireR;
222 int TLayerNo, TWireNo,TSignalLayerNo,fSignalLayer[50];
223 int i,wireNo, firstWire,segmentNo,signalLayer;
224 double length, nomphi, phi, r,nomaadiv2, aadiv2,
225 nomaa, aa, nomrotateCell, rotateCell, wlength, innerR, outR, z;
226 std::string layer, name, line;
227
228 vector<int> WireNo, fSegmentNo;
229 vector<double> wLength, R, nomPhi, Phi,nomadiv2, adiv2, First,
230 noma, a, nomebusing, ebusing, nomsigma, sigma, nomdeltaz, deltaz,
231 nomRotate, Rotate, fLength, fInnerR, fOutR, fZ;
232 vector<string> LayerName, fName;
233 vector<double> Sx,Sy,Sz,Rx,Ry,Rz;
234 double tmp1,tmp2,tmp3,tmp4,tmp5,tmp6;
235
236 /// read alignment file
237
238 if (m_readAlignParDataBase && m_updataalign) {
239 ReadAliParDataBase(Sx, Sy, Sz, Rx, Ry, Rz);
240 } else {
241 cout << "Read alignment parameters from file:" << alignFilePath.c_str() << endl;
242 getline(alignFile, line);
243 alignFile.seekg(1,ios::cur);
244 for(int k=0;k<16;k++){
245 alignFile>>line>>tmp1>>tmp2>>tmp3>>tmp4>>tmp5>>tmp6;
246 //changed by luot 10.03.12
247 //if (!m_updataalign) {
248 // tmp1 = 0.0;
249 // tmp2 = 0.0;
250 // tmp3 = 0.0;
251 // tmp4 = 0.0;
252 // tmp5 = 0.0;
253 // tmp6 = 0.0;
254 //}
255 Sx.push_back(m_wholeShiftX+tmp1);
256 Sy.push_back(m_wholeShiftY+tmp2);
257 Sz.push_back(m_wholeShiftZ+tmp3);
258 Rx.push_back(m_wholeRotatX+tmp4);
259 Ry.push_back(m_wholeRotatY+tmp5);
260 Rz.push_back(m_wholeRotatZ+tmp6);
261 }
262 /// test the service
263 /// for(int m=0;m<16;m++) std::cout<<" Rz["<<m<<"] is: "<<Rz[m]<<std::endl;
264 }
265
266 getline(inFile, line);
267 inFile>>TLayerNo>>TWireNo>>TSignalLayerNo>>signalWireR>>fieldWireR;
268 inFile.seekg(1,ios::cur);
269 getline(inFile, line);
270 for( i=0; i<TSignalLayerNo; i++){
271 inFile>>signalLayer;
272 /// signalLayer is the index of this signalLayer in all the layers
273 fSignalLayer[i]=signalLayer-1;
274 }
275
276 inFile.seekg(1,ios::cur);
277 getline(inFile, line);
278 getline(inFile, line);
279
280 int i_east,i_west;
281 double delta_rotateCell;
282 /// in this loop we read the file and store them into according vectors
283 for( i=0; i<TLayerNo; i++){
284 i_east=getAlignParIndexEast(i);
285 i_west=getAlignParIndexWest(i);
286
287 inFile>>layer>>wireNo>>wlength>>r>>nomphi>>firstWire>>nomrotateCell;
288 getline(inFile, line);
289
290 /// wireNo is the total wire number of this layer, including both signal
291 /// wire& field wire
292 nomaadiv2=2*pi*nomrotateCell/wireNo;
293 //aadiv2=2*pi*rotateCell/wireNo+0.5*(Rz[i_west]-Rz[i_east]);
294 nomaa=2*r*sin(nomaadiv2);
295 delta_rotateCell=(Rz[i_west]-Rz[i_east])*wireNo/(4*pi);
296 rotateCell=nomrotateCell+delta_rotateCell;
297 aadiv2=2*pi*rotateCell/wireNo;
298 aa=2*r*sin(aadiv2);
299
300 double shiftPhi=twopi/wireNo;
301 nomphi*=deg;
302 if(nomphi<0) nomphi += shiftPhi; //By definition, phi of first wire must >=0
303 phi=nomphi+Rz[i_east];
304 LayerName.push_back(layer);
305 WireNo.push_back(wireNo);
306 wLength.push_back(wlength);
307 R.push_back(r);
308 nomPhi.push_back(nomphi);
309 Phi.push_back(phi);
310 //nomadiv2.push_back(nomaadiv2);
311 //adiv2.push_back(aadiv2);
312 First.push_back(firstWire);
313 //noma.push_back(nomaa);
314 a.push_back(aa);
315 /// CLHEP says that:
316 /// static const double radian = 1.;
317 /// static const double degree = (3.14159265358979323846/180.0)*radian
318 /// static const double deg = degree;
319
320 nomebusing.push_back(atan(nomaa/wlength));
321 ebusing.push_back(atan(aa/wlength));
322
323 //if(aa==0.0) sigma.push_back(0.0);
324 //else sigma.push_back(0.13*wlength/aa);
325 //nomdeltaz.push_back(r*(1.0-cos(nomaadiv2)));
326 //deltaz.push_back(r*(1.0-cos(aadiv2)));
327 nomRotate.push_back(nomrotateCell);
328 Rotate.push_back(rotateCell);
329 }
330
331 getline(inFile, line);
332 inFile>>segmentNo;
333 inFile.seekg(1,ios::cur);
334 getline(inFile, line);
335 getline(inFile, line);
336
337 for(i=0; i<segmentNo; i++){
338 inFile>>length>>innerR>>outR>>z>>name;
339 getline(inFile,line);
340 fLength.push_back(length);
341 fInnerR.push_back(innerR);
342 fOutR.push_back(outR);
343 fZ.push_back(z);
344 fName.push_back(name);
345 }
346
347 /// Variables used here are as follows,
348 /// MdcGeoGeneral hlh;
349 /// MdcGeoSuper* suph;
350 /// MdcGeoLayer* lh;
351 /// MdcGeoWire* wh;
352 /// MdcGeoEnd* end;
353
354 //MdcGeoGeneral include signal wire and field wire
355 MdcGeoGeneral hlh;
356 double a1,a2,a3,nomphi0,phi0,cellphi;
357 double noma1,noma2,noma3;
358 for(i=0;i<TLayerNo;i++){
359 i_east=getAlignParIndexEast(i);
360 i_west=getAlignParIndexWest(i);
361 hlh.Id(i);
362 hlh.LayerName(LayerName[i]);
363 hlh.Radius(R[i]);
364 hlh.Length(wLength[i]);
365 hlh.NCell(WireNo[i]/2);
366 //east endcap rotates away from designed position by Phi considering misalignment,
367 //nomPhi without considering misalignment
368 hlh.Offset(Phi[i]);
369 hlh.nomOffset(nomPhi[i]);
370 hlh.Phi(Phi[i]);
371 hlh.nomPhi(nomPhi[i]);
372 hlh.First((int)First[i]);
373 HepPoint3D offF(Sx[i_west],Sy[i_west],Sz[i_west]);
374 HepPoint3D offB(Sx[i_east],Sy[i_east],Sz[i_east]);
375 hlh.OffF(offF);
376 hlh.OffB(offB);
377 hlh.TwistF(Rz[i_west]);
378 hlh.TwistB(Rz[i_east]);
379 //west endcap rotates away from designed position by Rotate considering misalignment,
380 //nomRotate without considering misalignment
381 hlh.Shift(Rotate[i]);
382 hlh.nomShift(nomRotate[i]);
383 hlh.SxEast(Sx[i_east]);
384 hlh.SxWest(Sx[i_west]);
385 hlh.SyEast(Sy[i_east]);
386 hlh.SyWest(Sy[i_west]);
387 hlh.SzEast(Sz[i_east]);
388 hlh.SzWest(Sz[i_west]);
389
390 hlh.RxEast(Rx[i_east]);
391 hlh.RxWest(Rx[i_west]);
392 hlh.RyEast(Ry[i_east]);
393 hlh.RyWest(Ry[i_west]);
394 hlh.RzEast(Rz[i_east]);
395 hlh.RzWest(Rz[i_west]);
396 fGenerals.push_back(hlh);
397 }
398
399 MdcGeoSuper *suph;
400 for(i=0;i<TSignalLayerNo;i++){
401 //MdcGeoLayer only represent signal wires
402 MdcGeoLayer *lh = new MdcGeoLayer();
403 // lh is a signalwire layer
404 int lyr=fSignalLayer[i];
405 // the value of lyr starts from 1,because the first layer of MDC is field
406 // wire
407 i_east=getAlignParIndexEast(lyr);
408 i_west=getAlignParIndexWest(lyr);
409 lh->SLayer(lyr+1);
410 if(i%4==0){
411 suph = new MdcGeoSuper();
412 suph->LowerR(R[lyr-1]);
413 if(i==40){
414 suph->UpperR(R[lyr+5]);
415 }
416 else {
417 suph->UpperR(R[lyr+7]);
418 }
419 switch(i/4){
420 case 0: suph->Type(-1); break;
421 case 1: suph->Type( 1); break;
422 case 2: suph->Type( 0); break;
423 case 3: suph->Type( 0); break;
424 case 4: suph->Type( 0); break;
425 case 5: suph->Type(-1); break;
426 case 6: suph->Type( 1); break;
427 case 7: suph->Type(-1); break;
428 case 8: suph->Type( 1); break;
429 case 9: suph->Type( 0); break;
430 case 10: suph->Type(0); break;
431 default: break;
432 }
433 suph->Id(i/4);
434 fSupers.push_back(suph);
435 }
436 /// attention, it is 2*twopi means 4*pi
437 cellphi=2*twopi/WireNo[lyr];
438 /// calculate the phi angle of the sense wire of the first cell of the layer
439 phi0=Phi[lyr]+abs(First[lyr]-1)*(cellphi/2);
440 nomphi0=nomPhi[lyr]+abs(First[lyr]-1)*(cellphi/2);
441 lh->Id(i);
442 int wircount=0;
443 for(int wk=0; wk<i; wk++){
444 int lyrwk=fSignalLayer[wk];
445 wircount=wircount + WireNo[lyrwk]/2;
446 };
447 /// Wirst is the the global number of the first signal wire of this layer
448 /// in all the signal wire
449 lh->Wirst(wircount);
450 lh->Slant(ebusing[lyr]);
451 lh->nomSlant(nomebusing[lyr]);
452 lh->Radius(R[lyr]);
453 lh->Length(wLength[lyr]);
454 lh->RCSiz1(R[lyr]-R[lyr-1]);
455 lh->RCSiz2(R[lyr+1]-R[lyr]);
456 lh->PCSiz(R[lyr]*cellphi);
457 lh->NCell(WireNo[lyr]/2);
458 lh->Offset(phi0);
459 lh->Shift(Rotate[lyr]);
460 lh->nomOffset(nomphi0);
461 lh->nomShift(nomRotate[lyr]);
462
463 vector<MdcGeoSuper*>::iterator it1 = fSupers.end()-1;
464 lh->Sup(* it1);
465 fLayers.push_back(lh);
466 int geo = -1;
467 if (i < 8) {
468 geo = i * 2 + 1;
469 } else if (i < 20) {
470 geo = i * 2 + 2;
471 } else if (i < 36) {
472 geo = i * 2 + 3;
473 } else {
474 geo = i * 2 + 4;
475 }
476 lh->Gen(geo);
477
478 for(int j=0;j<WireNo[lyr]/2;j++){
479 /// wh is signal wire
480 MdcGeoWire *wh = new MdcGeoWire();
481 const int wireId = wircount+j;
482 wh->Id(wireId);
483
484 /// caluculate misaligned wire position
485 a1 = R[lyr]*cos(phi0+j*cellphi)+Sx[i_east]+wirePosVec[wireId][0];
486 a2 = R[lyr]*sin(phi0+j*cellphi)+Sy[i_east]+wirePosVec[wireId][1];
487 a3 = wLength[lyr]/2+wirePosVec[wireId][2]+Sz[i_east];
488 //a3 = wLength[lyr]/2+wirePosVec[wireId][2];
489 //MDC rotate around x firstly;
490 double xb = a1;
491 double yb = a2 * cos(Rx[i_east]) + a3 * sin(Rx[i_east]);
492 double zb = a3 * cos(Rx[i_east]) - a2 * sin(Rx[i_east]);
493
494 //then MDC rotate around y;
495 double xbnew = xb * cos(Ry[i_east]) - zb * sin(Ry[i_east]);
496 double ybnew = yb;
497 double zbnew = xb * sin(Ry[i_east]) + zb * cos(Ry[i_east]);
498 ///
499 noma1 = R[lyr]*cos(nomphi0+j*cellphi);
500 noma2 = R[lyr]*sin(nomphi0+j*cellphi);
501 /// attention wLength[lyr] not the true length of the wire
502 noma3 = wLength[lyr]/2;
503 //HepPoint3D pb(a1,a2,a3);
504 HepPoint3D pb(xbnew,ybnew,zbnew);
505 HepPoint3D nompb(noma1,noma2,noma3);
506 HepPoint3D wireposb(wirePosVec[wireId][0],wirePosVec[wireId][1],wirePosVec[wireId][2]);
507 wh->Backward(pb);
508 wh->nomBackward(nompb);
509 wh->BWirePos(wireposb);
510
511 //double g1 = noma1*cos(Rz[i_east]) - noma2*sin(Rz[i_east]) + Sx[i_east] + wirePosVec[wireId][0];
512 //double g2 = noma2*cos(Rz[i_east]) + noma1*sin(Rz[i_east]) + Sy[i_east] + wirePosVec[wireId][1];
513 //double g3 = noma3 + wirePosVec[wireId][2];
514 //HepPoint3D backward_test(g1,g2,g3);
515
516 a1 = R[lyr]*cos(phi0+(j+lh->Shift())*cellphi)+Sx[i_west]+wirePosVec[wireId][3];
517 a2 = R[lyr]*sin(phi0+(j+lh->Shift())*cellphi)+Sy[i_west]+wirePosVec[wireId][4];
518 //a3 = -wLength[lyr]/2+wirePosVec[wireId][5];
519 a3 = -wLength[lyr]/2+wirePosVec[wireId][5]+Sz[i_west];
520 //MDC rotate around x firstly;
521 double xf = a1;
522 double yf = a2 * cos(Rx[i_west]) + a3 * sin(Rx[i_west]);
523 double zf = a3 * cos(Rx[i_west]) - a2 * sin(Rx[i_west]);
524
525 //then MDC rotate around y;
526 double xfnew = xf * cos(Ry[i_west]) - zf * sin(Ry[i_west]);
527 double yfnew = yf;
528 double zfnew = xf * sin(Ry[i_west]) + zf * cos(Ry[i_west]);
529
530 noma1 = R[lyr]*cos(nomphi0+(j+lh->nomShift())*cellphi);
531 noma2 = R[lyr]*sin(nomphi0+(j+lh->nomShift())*cellphi);
532 noma3 = -wLength[lyr]/2;
533 //HepPoint3D pf(a1,a2,a3);
534 HepPoint3D pf(xfnew, yfnew, zfnew);
535 HepPoint3D nompf(noma1,noma2,noma3);
536 HepPoint3D wireposf(wirePosVec[wireId][3],wirePosVec[wireId][4],wirePosVec[wireId][5]);
537 wh->Forward(pf);
538 wh->nomForward(nompf);
539 wh->FWirePos(wireposf);
540
541 //double f1 = noma1*cos(Rz[i_west]) - noma2*sin(Rz[i_west]) + Sx[i_west] + wirePosVec[wireId][3];
542 //double f2 = noma2*cos(Rz[i_west]) + noma1*sin(Rz[i_west]) + Sy[i_west] + wirePosVec[wireId][4];
543 //double f3 = noma3 + wirePosVec[wireId][5];
544 //HepPoint3D forward_test(f1,f2,f3);
545 /// test
546 /// std::cout<<" wireId: "<<wireId<<" backward position:"<< wh->Backward()<<" forward position: "<<wh->Forward()<<std::endl;
547 wh->Layer(i);
548 wh->Cell(j);
549 wh->Stat(0);
550 wh->Slant(ebusing[lyr]);
551 wh->nomSlant(nomebusing[lyr]);
552
553 /// set the tension and wire sag
554 wh->Tension(wireTensionVec[wireId]);
555 //std::cout<<" wh->Tension: "<< wh->Tension()<<std::endl;
556 //std::cout<<" wire sag: "<<wh->Sag(wireId)<<std::endl;
557
558 /// because of just before have fLayers.push_back(lh);
559 /// so use it2 = fLayers.end()-1 to get the MdcGeoLayer
560 /// which have just push_back into the stack
561 vector<MdcGeoLayer*>::iterator it2 = fLayers.end()-1;
562 wh->Lyr(* it2);
563 fWires.push_back(wh);
564 }
565 }
566
567 fMisc.OuterR(fOutR[1]);
568 fMisc.InnerR(fInnerR[2]);
569 fMisc.OuterTk(fOutR[1]-fInnerR[1]);
570 fMisc.InnerTk(fOutR[2]-fInnerR[2]);
571 fMisc.NSWire(fWires.size());
572 fMisc.NFWire(TWireNo-fWires.size());
573
574 fMisc.LayerNo(TLayerNo);
575 fMisc.WireNo(TWireNo);
576 fMisc.SLayerNo(TSignalLayerNo);
577 fMisc.SWireR(signalWireR);
578 fMisc.FWireR(fieldWireR);
579
580 for(i=0;i<segmentNo;i++){
581 MdcGeoEnd * end=new MdcGeoEnd();
582 end->Id(i);
583 end->Length(fLength[i]);
584 end->InnerR(fInnerR[i]);
585 end->OutR(fOutR[i]);
586 end->Z(fZ[i]);
587 end->Name(fName[i]);
588 fEnd.push_back(end);
589 }
590}
591
592void MdcGeomSvc::ReadTensionDataBase(std::vector<double> & wireTensionVec){
593 std::string fullPath = "/Calib/MdcAlign";
594 MsgStream log(messageService(), name());
595 log << MSG::INFO<<"ReadTensionDataBase() wireTensionPath = "<<fullPath<< endreq;
596 cout << "Read wireTension from Calibration Database!" << endl;
597
598 SmartDataPtr<CalibData::MdcAlignData> tension(m_pCalibDataSvc, fullPath);
599 if( ! tension){
600 log << MSG::ERROR << "can not get MdcAlignConst via SmartPtr"
601 << endreq;
602 return;
603 }
604 for(int i=0;i<6796;i++){
605 double tens = tension->gettension(i);
606 wireTensionVec.push_back(tens);
607 }
608}
609void MdcGeomSvc::ReadWirePosDataBase(std::vector<vector<double> > & wirePosVec){
610 std::string fullPath = "/Calib/MdcAlign";
611 MsgStream log(messageService(), name());
612 log << MSG::INFO<<"ReadWirePosDataBase() wirePositionPath = "<<fullPath<< endreq;
613
614 cout << "Read wirePosition from Calibration Database!" << endl;
615 SmartDataPtr<CalibData::MdcAlignData> wirepos(m_pCalibDataSvc, fullPath);
616 if( ! wirepos){
617 log << MSG::ERROR << "can not get MdcAlignConst via SmartPtr"
618 << endreq;
619 return;
620 }
621 double dxe,dye,dze,dxw,dyw,dzw;
622 for(int j=0; j<6796; j++){
623
624 dxe = wirepos->getdxWireEast(j);
625 dye = wirepos->getdyWireEast(j);
626 dze = wirepos->getdzWireEast(j);
627 dxw = wirepos->getdxWireWest(j);
628 dyw = wirepos->getdyWireWest(j);
629 dzw = wirepos->getdzWireWest(j);
630
631 wirePosVec[j].push_back(dxe);
632 wirePosVec[j].push_back(dye);
633 wirePosVec[j].push_back(dze);
634 wirePosVec[j].push_back(dxw);
635 wirePosVec[j].push_back(dyw);
636 wirePosVec[j].push_back(dzw);
637 }
638
639}
640void MdcGeomSvc::ReadAliParDataBase(vector<double>& Sx, vector<double>& Sy, vector<double>& Sz,
641 vector<double>& Rx, vector<double>& Ry, vector<double>& Rz){
642 MsgStream log(messageService(), name());
643 std::string fullPath = "/Calib/MdcAlign";
644 log << MSG::INFO<<"ReadAliParDataBase() alignParPath = "<<fullPath<< endreq;
645 cout << "Read alignment parameters from Calibration Database!" << endl;
646
647 SmartDataPtr<CalibData::MdcAlignData> alignpar(m_pCalibDataSvc, fullPath);
648 if( ! alignpar){
649 log << MSG::ERROR << "can not get MdcAlignConst via SmartPtr"
650 << endreq;
651 return;
652 }
653 double tmp1,tmp2,tmp3,tmp4,tmp5,tmp6;
654 for(int k=0;k<16;k++){
655 tmp1 = alignpar->getdxEP(k);
656 tmp2 = alignpar->getdyEP(k);
657 tmp3 = alignpar->getdzEP(k);
658 tmp4 = alignpar->getrxEP(k);
659 tmp5 = alignpar->getryEP(k);
660 tmp6 = alignpar->getrzEP(k);
661
662 Sx.push_back(m_wholeShiftX+tmp1);
663 Sy.push_back(m_wholeShiftY+tmp2);
664 Sz.push_back(m_wholeShiftZ+tmp3);
665 Rx.push_back(m_wholeRotatX+tmp4);
666 Ry.push_back(m_wholeRotatY+tmp5);
667 Rz.push_back(m_wholeRotatZ+tmp6);
668 }
669
670}
672{
673 return fWires.size();
674}
675
677{
678 return fLayers.size();
679}
680
682{
683 return fSupers.size();
684}
685
687{
688 return fGenerals.size();
689}
690
692{
693 return fEnd.size();
694}
695
697
698const int MdcGeomSvc::getAlignParIndexEast(int lyr) const{
699 /// the Mdc small endplate east
700 if((lyr>=0) && (lyr<=16)) return 0;
701 /// the 1th ring east
702 else if((lyr>=17) && (lyr<=20)) return 1;
703 /// the 2th ring east
704 else if((lyr>=21) && (lyr<=24)) return 2;
705 /// the 3th ring east
706 else if((lyr>=25) && (lyr<=28)) return 3;
707 /// the 4th ring east
708 else if((lyr>=29) && (lyr<=32)) return 4;
709 /// the 5th ring east
710 else if((lyr>=33) && (lyr<=36)) return 5;
711 /// the 6th ring east
712 else if((lyr>=37) && (lyr<=41)) return 6;
713 /// the Mdc big endplate east
714 else if(lyr>=42) return 7;
715 else std::cout<<" Hi, ERROR OCCUR !!!"<<std::endl;
716 return -1;
717}
718
719
720const int MdcGeomSvc::getAlignParIndexWest(int lyr) const{
721 /// the Mdc small endplate west
722 if((lyr>=0) && (lyr<=16)) return 8;
723 /// the 1th ring west
724 else if((lyr>=17) && (lyr<=20)) return 9;
725 /// the 2th ring west
726 else if((lyr>=21) && (lyr<=24)) return 10;
727 /// the 3th ring west
728 else if((lyr>=25) && (lyr<=28)) return 11;
729 /// the 4th ring west
730 else if((lyr>=29) && (lyr<=32)) return 12;
731 /// the 5th ring west
732 else if((lyr>=33) && (lyr<=36)) return 13;
733 /// the 6th ring west
734 else if((lyr>=37) && (lyr<=41)) return 14;
735 /// the Mdc big endplate west
736 else if(lyr>=42) return 15;
737 else std::cout<<" Hi, ERROR OCCUR !!!"<<std::endl;
738 return -1;
739}
740
741
742/// this handle function is prepared for special use
743void MdcGeomSvc::handle(const Incident& inc){
744 MsgStream log( messageService(), name() );
745 log << MSG::DEBUG << "handle: " << inc.type() << endreq;
746 IDataProviderSvc* m_eventSvc;
747 Gaudi::svcLocator()->service("EventDataSvc", m_eventSvc, true);
748 SmartDataPtr<Event::EventHeader> eventHeader(m_eventSvc,"/Event/EventHeader");
749 if (!eventHeader) {
750 log << MSG::FATAL << "Could not find Event Header" << endreq;
751 }
752 if (m_updataalign) return;
753 if (inc.type() == "NewRun" ){
754 log << MSG::DEBUG << "Begin Event" << endreq;
755 clean();
756 m_updataalign = true;
757 if(m_nomcalignment&&m_mindex==0) {
758 int RunNo=eventHeader->runNumber();
759 if(RunNo<0) m_readAlignParDataBase=false ;
760 else m_readAlignParDataBase=true;
761 m_mindex+=1;
762 cout<<"m__RunNo="<<RunNo<<"m_mindex="<<m_mindex<<endl;
763 }
764 //std::cout<<"############"<<m_readAlignParDataBase<<std::endl;
765 ReadFilePar();
766 }
767}
768
769const MdcGeoWire * const
770 MdcGeomSvc::Wire(unsigned id){
771 if (id < fWires.size())
772 return fWires[id];
773
774 return 0;
775 }
776
777const MdcGeoWire * const
778 MdcGeomSvc::Wire(unsigned lyrid, unsigned wirid){
779 if ((lyrid <fLayers.size()) && ((int) wirid < Layer(lyrid)->NCell()))
780 return fWires[Layer(lyrid)->Wirst() + wirid];
781
782 return 0;
783 }
784
785const MdcGeoLayer * const
786 MdcGeomSvc::Layer(unsigned id){
787 if (id < fLayers.size())
788 return fLayers[id];
789
790 return 0;
791 }
792
793const MdcGeoSuper * const
795 if (id < fSupers.size())
796 return fSupers[id];
797
798 return 0;
799 }
800
801const MdcGeoGeneral * const
803 if (id < fGenerals.size())
804 return & fGenerals[id];
805
806 return 0;
807 }
808
809const MdcGeoMisc * const
811 return & fMisc;
812}
813
814const MdcGeoEnd * const
815 MdcGeomSvc::End(unsigned id){
816 if (id < fEnd.size())
817 return fEnd[id];
818
819 return 0;
820 }
821
822
824
825 return m_doSag;
826}
827
828
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
double Phi(RecMdcKalTrack *trk)
double length
const int wireNo
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
#define DBL_MAX
Definition: KalFitAlg.h:13
double Length(void) const
Definition: MdcGeoEnd.h:17
string Name(void) const
Definition: MdcGeoEnd.h:21
double InnerR(void) const
Definition: MdcGeoEnd.h:18
double OutR(void) const
Definition: MdcGeoEnd.h:19
double Z(void) const
Definition: MdcGeoEnd.h:20
int Id(void) const
Definition: MdcGeoEnd.h:16
double Length(void) const
int Id(void) const
double nomShift(void) const
double TwistB(void) const
double SzWest(void) const
HepPoint3D OffB(void) const
int NCell(void) const
double nomPhi(void) const
double Phi(void) const
double RxEast(void) const
double SxEast(void) const
double TwistF(void) const
double Offset(void) const
double Shift(void) const
double SyWest(void) const
double SxWest(void) const
string LayerName(void) const
double RxWest(void) const
double nomOffset(void) const
double RyWest(void) const
HepPoint3D OffF(void) const
int First(void) const
double SzEast(void) const
double RzWest(void) const
double RzEast(void) const
double SyEast(void) const
double RyEast(void) const
double Radius(void) const
int Id(void) const
Definition: MdcGeoLayer.h:155
double Radius(void) const
Definition: MdcGeoLayer.h:160
double Slant(void) const
Definition: MdcGeoLayer.h:158
int Gen(void) const
Definition: MdcGeoLayer.h:175
double RCSiz2(void) const
Definition: MdcGeoLayer.h:163
double Length(void) const
Definition: MdcGeoLayer.h:161
double PCSiz(void) const
Definition: MdcGeoLayer.h:164
MdcGeoSuper * Sup(void) const
Definition: MdcGeoLayer.h:174
int SLayer(void) const
Definition: MdcGeoLayer.h:156
double nomShift(void) const
Definition: MdcGeoLayer.h:169
double Shift(void) const
Definition: MdcGeoLayer.h:167
double nomSlant(void) const
Definition: MdcGeoLayer.h:159
double RCSiz1(void) const
Definition: MdcGeoLayer.h:162
int NCell(void) const
Definition: MdcGeoLayer.h:165
int Wirst(void) const
Definition: MdcGeoLayer.h:157
double Offset(void) const
Definition: MdcGeoLayer.h:166
double nomOffset(void) const
Definition: MdcGeoLayer.h:168
double FWireR(void) const
Definition: MdcGeoMisc.h:53
int NSWire(void) const
Definition: MdcGeoMisc.h:46
double OuterTk(void) const
Definition: MdcGeoMisc.h:44
int LayerNo(void) const
Definition: MdcGeoMisc.h:49
double InnerR(void) const
Definition: MdcGeoMisc.h:43
int WireNo(void) const
Definition: MdcGeoMisc.h:50
double OuterR(void) const
Definition: MdcGeoMisc.h:42
int NFWire(void) const
Definition: MdcGeoMisc.h:47
double SWireR(void) const
Definition: MdcGeoMisc.h:52
double InnerTk(void) const
Definition: MdcGeoMisc.h:45
int SLayerNo(void) const
Definition: MdcGeoMisc.h:51
double LowerR(void) const
Definition: MdcGeoSuper.h:34
int Id(void) const
Definition: MdcGeoSuper.h:32
double UpperR(void) const
Definition: MdcGeoSuper.h:33
int Type(void) const
Definition: MdcGeoSuper.h:35
HepPoint3D FWirePos(void) const
Definition: MdcGeoWire.h:131
double Slant(void) const
Definition: MdcGeoWire.h:132
int Cell(void) const
Definition: MdcGeoWire.h:137
HepPoint3D BWirePos(void) const
Definition: MdcGeoWire.h:130
int Stat(void) const
Definition: MdcGeoWire.h:139
MdcGeoLayer * Lyr(void) const
Definition: MdcGeoWire.h:140
HepPoint3D nomForward(void) const
Definition: MdcGeoWire.h:134
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
int Id(void) const
Definition: MdcGeoWire.h:127
HepPoint3D nomBackward(void) const
Definition: MdcGeoWire.h:133
int Layer(void) const
Definition: MdcGeoWire.h:138
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128
double Tension(void) const
Definition: MdcGeoWire.h:136
double nomSlant(void) const
Definition: MdcGeoWire.h:135
static bool m_readAlignParDataBase
Definition: MdcGeomSvc.h:55
MdcGeomSvc(const std::string &name, ISvcLocator *svcloc)
Definition: MdcGeomSvc.cxx:30
void Dump()
Definition: MdcGeomSvc.cxx:696
virtual StatusCode initialize()
Definition: MdcGeomSvc.cxx:73
static bool getSagFlag(void)
Definition: MdcGeomSvc.cxx:823
const MdcGeoSuper *const SuperLayer(unsigned id)
Definition: MdcGeomSvc.cxx:794
const MdcGeoWire *const Wire(unsigned id)
Definition: MdcGeomSvc.cxx:770
const MdcGeoGeneral *const GeneralLayer(unsigned id)
Definition: MdcGeomSvc.cxx:802
void handle(const Incident &inc)
this handle function is prepared for special use
Definition: MdcGeomSvc.cxx:743
const int getSuperLayerSize()
Definition: MdcGeomSvc.cxx:681
virtual StatusCode finalize()
Definition: MdcGeomSvc.cxx:107
virtual StatusCode queryInterface(const InterfaceID &riid, void **ppvUnknown)
Definition: MdcGeomSvc.cxx:63
const MdcGeoEnd *const End(unsigned id)
Definition: MdcGeomSvc.cxx:815
const MdcGeoLayer *const Layer(unsigned id)
Definition: MdcGeomSvc.cxx:786
static bool m_doSag
Definition: MdcGeomSvc.h:54
const int getSegmentNo()
Definition: MdcGeomSvc.cxx:691
const int getWireSize()
Definition: MdcGeomSvc.cxx:671
const int getGeneralLayerSize()
Definition: MdcGeomSvc.cxx:686
const int getLayerSize()
Definition: MdcGeomSvc.cxx:676
const MdcGeoMisc *const Misc(void)
Definition: MdcGeomSvc.cxx:810
static bool m_nomcalignment
Definition: MdcGeomSvc.h:56
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition: TUtil.h:27
const float pi
Definition: vector3.h:133