CGEM BOSS 6.6.5.f
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 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
113MdcGeomSvc::~MdcGeomSvc(){
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}
671const int MdcGeomSvc::getWireSize()
672{
673 return fWires.size();
674}
675
676const int MdcGeomSvc::getLayerSize()
677{
678 return fLayers.size();
679}
680
681const int MdcGeomSvc::getSuperLayerSize()
682{
683 return fSupers.size();
684}
685
686const int MdcGeomSvc::getGeneralLayerSize()
687{
688 return fGenerals.size();
689}
690
691const int MdcGeomSvc::getSegmentNo()
692{
693 return fEnd.size();
694}
695
696void MdcGeomSvc::Dump(){}
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
794 MdcGeomSvc::SuperLayer(unsigned id){
795 if (id < fSupers.size())
796 return fSupers[id];
797
798 return 0;
799 }
800
801const MdcGeoGeneral * const
802 MdcGeomSvc::GeneralLayer(unsigned id){
803 if (id < fGenerals.size())
804 return & fGenerals[id];
805
806 return 0;
807 }
808
809const MdcGeoMisc * const
810MdcGeomSvc::Misc(void){
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
823bool MdcGeomSvc::getSagFlag(void) {
824
825 return m_doSag;
826}
827
828
Phi
double Phi(RecMdcKalTrack *trk)
Definition: CalibEventSelect.cxx:99
wireNo
const int wireNo
Definition: DedxCalibWireGain.cxx:18
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
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:696
MdcGeomSvc::initialize
virtual StatusCode initialize()
Definition: MdcGeomSvc.cxx:73
MdcGeomSvc::getSagFlag
static bool getSagFlag(void)
Definition: MdcGeomSvc.cxx:823
MdcGeomSvc::SuperLayer
const MdcGeoSuper *const SuperLayer(unsigned id)
Definition: MdcGeomSvc.cxx:794
MdcGeomSvc::Wire
const MdcGeoWire *const Wire(unsigned id)
Definition: MdcGeomSvc.cxx:770
MdcGeomSvc::GeneralLayer
const MdcGeoGeneral *const GeneralLayer(unsigned id)
Definition: MdcGeomSvc.cxx:802
MdcGeomSvc::handle
void handle(const Incident &inc)
this handle function is prepared for special use
Definition: MdcGeomSvc.cxx:743
MdcGeomSvc::getSuperLayerSize
const int getSuperLayerSize()
Definition: MdcGeomSvc.cxx:681
MdcGeomSvc::finalize
virtual StatusCode finalize()
Definition: MdcGeomSvc.cxx:107
MdcGeomSvc::queryInterface
virtual StatusCode queryInterface(const InterfaceID &riid, void **ppvUnknown)
Definition: MdcGeomSvc.cxx:63
MdcGeomSvc::End
const MdcGeoEnd *const End(unsigned id)
Definition: MdcGeomSvc.cxx:815
MdcGeomSvc::Layer
const MdcGeoLayer *const Layer(unsigned id)
Definition: MdcGeomSvc.cxx:786
MdcGeomSvc::getSegmentNo
const int getSegmentNo()
Definition: MdcGeomSvc.cxx:691
MdcGeomSvc::getWireSize
const int getWireSize()
Definition: MdcGeomSvc.cxx:671
MdcGeomSvc::getGeneralLayerSize
const int getGeneralLayerSize()
Definition: MdcGeomSvc.cxx:686
MdcGeomSvc::getLayerSize
const int getLayerSize()
Definition: MdcGeomSvc.cxx:676
MdcGeomSvc::Misc
const MdcGeoMisc *const Misc(void)
Definition: MdcGeomSvc.cxx:810
rb::R
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition: TUtil.h:27