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