CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcGeomSvc-00-01-42/src/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#include <math.h>
23
24
25
26// initialize static members
27bool MdcGeomSvc::m_doSag = true;
28bool MdcGeomSvc::m_readAlignParDataBase = true;
29bool MdcGeomSvc::m_nomcalignment=true;
30
31MdcGeomSvc::MdcGeomSvc( const std::string& name, ISvcLocator* svcloc ) : Service(name, svcloc) {
32 if(getenv("MDCGEOMSVCROOT")){
33 m_alignFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/MdcAlignPar.dat");
34 //std::cout<<" the MDC alignment file: "<<m_alignFilePath<<std::endl;
35
36 m_wirePosFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/WirePosCalib.dat");
37 //std::cout<<" the MDC wire position file: "<<m_wirePosFilePath<<std::endl;
38
39 m_wireTensionFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/mdcWireTension.dat");
40 //std::cout<<" the MDC wire tension file: "<<m_wireTensionFilePath<<std::endl;
41
42 }
43 else {
44 std::cout<<"A fatal error, contact wangjk..."<<std::endl;
45 }
46
47 declareProperty("doSag", m_doSag = true);
48 declareProperty("readAlignParDataBase", m_readAlignParDataBase = true);
49 declareProperty("mcnoalignment",m_nomcalignment=true);
50 declareProperty("wholeShiftX",m_wholeShiftX = 0.);
51 declareProperty("wholeShiftY",m_wholeShiftY = 0.);
52 declareProperty("wholeShiftZ",m_wholeShiftZ = 0.);
53 declareProperty("wholeRotatX",m_wholeRotatX = 0.);
54 declareProperty("wholeRotatY",m_wholeRotatY = 0.);
55 declareProperty("wholeRotatZ",m_wholeRotatZ = 0.);
56 declareProperty("alignFilePath",m_alignFilePath);
57 declareProperty("wirePosFilePath",m_wirePosFilePath);
58 declareProperty("wireTensionFilePath",m_wireTensionFilePath);
59 declareProperty("useCGEM",m_useCgem = true);
60 declareProperty("NeighborRange",m_NeighborRange = 1.);
61 declareProperty("testAdjacentIntersection",m_conf_testAdjacentIntersection = false);
62
63
64}
65
66StatusCode MdcGeomSvc::queryInterface (const InterfaceID& riid, void** ppvInterface ){
67
68 if ( IID_IMdcGeomSvc.versionMatch(riid) ) {
69 *ppvInterface = static_cast<IMdcGeomSvc*> (this);
70 } else {
71 return Service::queryInterface(riid, ppvInterface) ;
72 }
73 return StatusCode::SUCCESS;
74}
75
76StatusCode MdcGeomSvc::initialize ( ) {
77 MsgStream log(messageService(), name());
78 log << MSG::INFO << name() << ": Start of run initialisation" << endreq;
79
80 StatusCode sc = Service::initialize();
81 if ( sc.isFailure() ) return sc;
82 m_mindex=0;
83 m_updataalign = false;
84 IIncidentSvc* incsvc;
85 sc = service("IncidentSvc", incsvc);
86 int priority = 100;
87 if( sc.isSuccess() ){
88 incsvc -> addListener(this, "NewRun", priority);
89 }
90
91 // ReadFilePar(); // get geometry data from file SimUtil/dat/Mdc.txt
92 // Fill(); // get geometry data from Database
93
94 sc = service("CalibDataSvc", m_pCalibDataSvc, true);
95
96 if ( !sc.isSuccess() ) {
97 log << MSG::ERROR
98 << "Could not get IDataProviderSvc interface of CalibXmlCnvSvc"
99 << endreq;
100 return sc;
101 } else {
102 log << MSG::DEBUG
103 << "Retrieved IDataProviderSvc interface of CalibXmlCnvSvc"
104 << endreq;
105 }
106 ReadFilePar();
107 AddWireNeighbors();
108 return StatusCode::SUCCESS;
109}
110
111StatusCode MdcGeomSvc::finalize ( ) {
112 if (m_conf_testAdjacentIntersection)\
113 testAdjacentIntersection();
114 MsgStream log(messageService(), name());
115 log << MSG::INFO << name() << ": End of Run" << endreq;
116 return StatusCode::SUCCESS;
117}
118
119MdcGeomSvc::~MdcGeomSvc(){
120 for(vector<MdcGeoLayer*>::iterator it1 = fLayers.begin(); it1 != fLayers.end(); it1++) delete *it1;
121 for(vector<MdcGeoSuper*>::iterator it2 = fSupers.begin(); it2 != fSupers.end(); it2++) delete *it2;
122 for(vector<MdcGeoWire*>::iterator it3 = fWires.begin(); it3 != fWires.end(); it3++) delete *it3;
123 for(vector<MdcGeoEnd*>::iterator it4 = fEnd.begin(); it4 != fEnd.end(); it4++) delete *it4;
124 fGenerals.clear();
125 fWires.clear();
126 fLayers.clear();
127 fSupers.clear();
128 fEnd.clear();
129}
130
131void MdcGeomSvc::clean(){
132 for(vector<MdcGeoLayer*>::iterator it1 = fLayers.begin(); it1 != fLayers.end(); it1++) delete *it1;
133 for(vector<MdcGeoSuper*>::iterator it2 = fSupers.begin(); it2 != fSupers.end(); it2++) delete *it2;
134 for(vector<MdcGeoWire*>::iterator it3 = fWires.begin(); it3 != fWires.end(); it3++) delete *it3;
135 for(vector<MdcGeoEnd*>::iterator it4 = fEnd.begin(); it4 != fEnd.end(); it4++) delete *it4;
136 fGenerals.clear();
137 fWires.clear();
138 fLayers.clear();
139 fSupers.clear();
140 fEnd.clear();
141}
142
143
144void MdcGeomSvc::ReadFilePar(){
145 std::string geometryFilePath = getenv("MDCSIMROOT");
146 if(m_useCgem) geometryFilePath += "/dat/Mdc_cgem.txt";
147 else geometryFilePath += "/dat/Mdc.txt";
148 std::ifstream inFile(geometryFilePath.c_str() );
149
150 if(!inFile.good()){
151 std::cout<<"Error, mdc parameters file not exist"<<std::endl;
152 return;
153 }
154 std::string alignFilePath;
155 std::string wirePosFilePath;
156 std::string wireTensionFilePath;
157 if (!m_updataalign) {
158 alignFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/MdcAlignPar.dat");
159 wirePosFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/WirePosCalib.dat");
160 wireTensionFilePath = std::string(getenv("MDCGEOMSVCROOT"))+std::string("/share/mdcWireTension.dat");
161 } else {
162 alignFilePath = m_alignFilePath;
163 wirePosFilePath = m_wirePosFilePath;
164 wireTensionFilePath = m_wireTensionFilePath;
165 }
166 std::ifstream alignFile(alignFilePath.c_str());
167 std::ifstream wirePosFile(wirePosFilePath.c_str());
168 std::ifstream wireTensionFile(wireTensionFilePath.c_str());
169
170 if(!wirePosFile){
171 std::cout<<"Error, mdc wire position file not exist..."<<std::endl;
172 }
173
174 if(!wireTensionFile){
175 std::cout<<"Error, mdc wire tension file not exist..."<<std::endl;
176 }
177
178 std::vector<double> wireTensionVec;
179 //cout << __FILE__ << " " << __LINE__ << " readdb = " << m_readAlignParDataBase
180 // << " update align = " << m_updataalign << endl;
181 if (m_readAlignParDataBase && m_updataalign) {
182 ReadTensionDataBase(wireTensionVec);
183 }else {
184 cout << "Read wireTension from file:" << wireTensionFilePath.c_str() << endl;
185 int idd;
186 double tension;
187 for(int ii=0; ii<6796; ii++){
188 wireTensionFile>>idd>>tension;
189
190 if(getSagFlag()){
191 wireTensionVec.push_back(tension);
192 }
193 else {
194 tension = DBL_MAX;
195 wireTensionVec.push_back(tension);
196 }
197 }
198 }
199
200 std::vector<vector<double> > wirePosVec(6796);
201 if (m_readAlignParDataBase && m_updataalign) {
202 ReadWirePosDataBase(wirePosVec);
203 } else {
204 cout << "Read wirePosition from file:" << wirePosFilePath.c_str() << endl;
205 double dxe,dye,dze,dxw,dyw,dzw;
206 int wid;
207 std::string title;
208 getline(wirePosFile, title);
209 wirePosFile.seekg(1,ios::cur);
210 for(int j=0; j<6796; j++){
211 wirePosFile>>wid>>dxe>>dye>>dze>>dxw>>dyw>>dzw;
212 wirePosVec[j].push_back(dxe);
213 wirePosVec[j].push_back(dye);
214 wirePosVec[j].push_back(dze);
215 wirePosVec[j].push_back(dxw);
216 wirePosVec[j].push_back(dyw);
217 wirePosVec[j].push_back(dzw);
218 /// test
219 /*
220 for(int jj=0;jj<6;jj++){
221 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;
222 }
223 */
224 }
225 }
226
227 double signalWireR, fieldWireR;
228 int TLayerNo, TWireNo,TSignalLayerNo,fSignalLayer[50];
229 int i,wireNo, firstWire,segmentNo,signalLayer;
230 double length, nomphi, phi, r,nomaadiv2, aadiv2,
231 nomaa, aa, nomrotateCell, rotateCell, wlength, innerR, outR, z;
232 std::string layer, name, line;
233
234 vector<int> WireNo, fSegmentNo;
235 vector<double> wLength, R, nomPhi, Phi,nomadiv2, adiv2, First,
236 noma, a, nomebusing, ebusing, nomsigma, sigma, nomdeltaz, deltaz,
237 nomRotate, Rotate, fLength, fInnerR, fOutR, fZ;
238 vector<string> LayerName, fName;
239 vector<double> Sx,Sy,Sz,Rx,Ry,Rz;
240 double tmp1,tmp2,tmp3,tmp4,tmp5,tmp6;
241
242 /// read alignment file
243
244 if (m_readAlignParDataBase && m_updataalign) {
245 ReadAliParDataBase(Sx, Sy, Sz, Rx, Ry, Rz);
246 } else {
247 cout << "Read alignment parameters from file:" << alignFilePath.c_str() << endl;
248 getline(alignFile, line);
249 alignFile.seekg(1,ios::cur);
250 for(int k=0;k<16;k++){
251 alignFile>>line>>tmp1>>tmp2>>tmp3>>tmp4>>tmp5>>tmp6;
252 //changed by luot 10.03.12
253 //if (!m_updataalign) {
254 // tmp1 = 0.0;
255 // tmp2 = 0.0;
256 // tmp3 = 0.0;
257 // tmp4 = 0.0;
258 // tmp5 = 0.0;
259 // tmp6 = 0.0;
260 //}
261 Sx.push_back(m_wholeShiftX+tmp1);
262 Sy.push_back(m_wholeShiftY+tmp2);
263 Sz.push_back(m_wholeShiftZ+tmp3);
264 Rx.push_back(m_wholeRotatX+tmp4);
265 Ry.push_back(m_wholeRotatY+tmp5);
266 Rz.push_back(m_wholeRotatZ+tmp6);
267 }
268 /// test the service
269 /// for(int m=0;m<16;m++) std::cout<<" Rz["<<m<<"] is: "<<Rz[m]<<std::endl;
270 }
271
272 getline(inFile, line);
273 inFile>>TLayerNo>>TWireNo>>TSignalLayerNo>>signalWireR>>fieldWireR;
274 inFile.seekg(1,ios::cur);
275 getline(inFile, line);
276 for( i=0; i<TSignalLayerNo; i++){
277 inFile>>signalLayer;
278 /// signalLayer is the index of this signalLayer in all the layers
279 fSignalLayer[i]=signalLayer-1;
280 }
281
282 inFile.seekg(1,ios::cur);
283 getline(inFile, line);
284 getline(inFile, line);
285
286 int i_east,i_west;
287 double delta_rotateCell;
288 /// in this loop we read the file and store them into according vectors
289 for( i=0; i<TLayerNo; i++){
290 i_east=getAlignParIndexEast(i);
291 i_west=getAlignParIndexWest(i);
292
293 inFile>>layer>>wireNo>>wlength>>r>>nomphi>>firstWire>>nomrotateCell;
294 getline(inFile, line);
295
296 /// wireNo is the total wire number of this layer, including both signal
297 /// wire& field wire
298 nomaadiv2=2*pi*nomrotateCell/wireNo;
299 //aadiv2=2*pi*rotateCell/wireNo+0.5*(Rz[i_west]-Rz[i_east]);
300 nomaa=2*r*sin(nomaadiv2);
301 delta_rotateCell=(Rz[i_west]-Rz[i_east])*wireNo/(4*pi);
302 rotateCell=nomrotateCell+delta_rotateCell;
303 aadiv2=2*pi*rotateCell/wireNo;
304 aa=2*r*sin(aadiv2);
305
306 double shiftPhi=twopi/wireNo;
307 nomphi*=deg;
308 if(nomphi<0) nomphi += shiftPhi; //By definition, phi of first wire must >=0
309 phi=nomphi+Rz[i_east];
310 LayerName.push_back(layer);
311 WireNo.push_back(wireNo);
312 wLength.push_back(wlength);
313 R.push_back(r);
314 nomPhi.push_back(nomphi);
315 Phi.push_back(phi);
316 //nomadiv2.push_back(nomaadiv2);
317 //adiv2.push_back(aadiv2);
318 First.push_back(firstWire);
319 //noma.push_back(nomaa);
320 a.push_back(aa);
321 /// CLHEP says that:
322 /// static const double radian = 1.;
323 /// static const double degree = (3.14159265358979323846/180.0)*radian
324 /// static const double deg = degree;
325
326 nomebusing.push_back(atan(nomaa/wlength));
327 ebusing.push_back(atan(aa/wlength));
328
329 //if(aa==0.0) sigma.push_back(0.0);
330 //else sigma.push_back(0.13*wlength/aa);
331 //nomdeltaz.push_back(r*(1.0-cos(nomaadiv2)));
332 //deltaz.push_back(r*(1.0-cos(aadiv2)));
333 nomRotate.push_back(nomrotateCell);
334 Rotate.push_back(rotateCell);
335 }
336
337 getline(inFile, line);
338 inFile>>segmentNo;
339 inFile.seekg(1,ios::cur);
340 getline(inFile, line);
341 getline(inFile, line);
342
343 for(i=0; i<segmentNo; i++){
344 inFile>>length>>innerR>>outR>>z>>name;
345 getline(inFile,line);
346 fLength.push_back(length);
347 fInnerR.push_back(innerR);
348 fOutR.push_back(outR);
349 fZ.push_back(z);
350 fName.push_back(name);
351 }
352
353 /// Variables used here are as follows,
354 /// MdcGeoGeneral hlh;
355 /// MdcGeoSuper* suph;
356 /// MdcGeoLayer* lh;
357 /// MdcGeoWire* wh;
358 /// MdcGeoEnd* end;
359
360 //MdcGeoGeneral include signal wire and field wire
361 MdcGeoGeneral hlh;
362 double a1,a2,a3,nomphi0,phi0,cellphi;
363 double noma1,noma2,noma3;
364 for(i=0;i<TLayerNo;i++){
365 i_east=getAlignParIndexEast(i);
366 i_west=getAlignParIndexWest(i);
367 hlh.Id(i);
368 hlh.LayerName(LayerName[i]);
369 hlh.Radius(R[i]);
370 hlh.Length(wLength[i]);
371 hlh.NCell(WireNo[i]/2);
372 //east endcap rotates away from designed position by Phi considering misalignment,
373 //nomPhi without considering misalignment
374 hlh.Offset(Phi[i]);
375 hlh.nomOffset(nomPhi[i]);
376 hlh.Phi(Phi[i]);
377 hlh.nomPhi(nomPhi[i]);
378 hlh.First((int)First[i]);
379 HepPoint3D offF(Sx[i_west],Sy[i_west],Sz[i_west]);
380 HepPoint3D offB(Sx[i_east],Sy[i_east],Sz[i_east]);
381 hlh.OffF(offF);
382 hlh.OffB(offB);
383 hlh.TwistF(Rz[i_west]);
384 hlh.TwistB(Rz[i_east]);
385 //west endcap rotates away from designed position by Rotate considering misalignment,
386 //nomRotate without considering misalignment
387 hlh.Shift(Rotate[i]);
388 hlh.nomShift(nomRotate[i]);
389 hlh.SxEast(Sx[i_east]);
390 hlh.SxWest(Sx[i_west]);
391 hlh.SyEast(Sy[i_east]);
392 hlh.SyWest(Sy[i_west]);
393 hlh.SzEast(Sz[i_east]);
394 hlh.SzWest(Sz[i_west]);
395
396 hlh.RxEast(Rx[i_east]);
397 hlh.RxWest(Rx[i_west]);
398 hlh.RyEast(Ry[i_east]);
399 hlh.RyWest(Ry[i_west]);
400 hlh.RzEast(Rz[i_east]);
401 hlh.RzWest(Rz[i_west]);
402 fGenerals.push_back(hlh);
403 }
404
405 MdcGeoSuper *suph;
406 for(i=0;i<TSignalLayerNo;i++){
407 //MdcGeoLayer only represent signal wires
408 MdcGeoLayer *lh = new MdcGeoLayer();
409 // lh is a signalwire layer
410 int lyr=fSignalLayer[i];
411 // the value of lyr starts from 1,because the first layer of MDC is field
412 // wire
413 i_east=getAlignParIndexEast(lyr);
414 i_west=getAlignParIndexWest(lyr);
415 lh->SLayer(lyr+1);
416 if(i%4==0){
417 suph = new MdcGeoSuper();
418 suph->LowerR(R[lyr-1]);
419 if(i==40){
420 suph->UpperR(R[lyr+5]);
421 }
422 else {
423 suph->UpperR(R[lyr+7]);
424 }
425 switch(i/4){
426 case 0: suph->Type(-1); break;
427 case 1: suph->Type( 1); break;
428 case 2: suph->Type( 0); break;
429 case 3: suph->Type( 0); break;
430 case 4: suph->Type( 0); break;
431 case 5: suph->Type(-1); break;
432 case 6: suph->Type( 1); break;
433 case 7: suph->Type(-1); break;
434 case 8: suph->Type( 1); break;
435 case 9: suph->Type( 0); break;
436 case 10: suph->Type(0); break;
437 default: break;
438 }
439 suph->Id(i/4);
440 fSupers.push_back(suph);
441 }
442 /// attention, it is 2*twopi means 4*pi
443 cellphi=2*twopi/WireNo[lyr];
444 /// calculate the phi angle of the sense wire of the first cell of the layer
445 phi0=Phi[lyr]+abs(First[lyr]-1)*(cellphi/2);
446 nomphi0=nomPhi[lyr]+abs(First[lyr]-1)*(cellphi/2);
447 lh->Id(i);
448 int wircount=0;
449 for(int wk=0; wk<i; wk++){
450 int lyrwk=fSignalLayer[wk];
451 wircount=wircount + WireNo[lyrwk]/2;
452 };
453 /// Wirst is the the global number of the first signal wire of this layer
454 /// in all the signal wire
455 lh->Wirst(wircount);
456 lh->Slant(ebusing[lyr]);
457 lh->nomSlant(nomebusing[lyr]);
458 lh->Radius(R[lyr]);
459 lh->Length(wLength[lyr]);
460 lh->RCSiz1(R[lyr]-R[lyr-1]);
461 lh->RCSiz2(R[lyr+1]-R[lyr]);
462 lh->PCSiz(R[lyr]*cellphi);
463 lh->NCell(WireNo[lyr]/2);
464 lh->Offset(phi0);
465 lh->Shift(Rotate[lyr]);
466 lh->nomOffset(nomphi0);
467 lh->nomShift(nomRotate[lyr]);
468
469 vector<MdcGeoSuper*>::iterator it1 = fSupers.end()-1;
470 lh->Sup(* it1);
471 fLayers.push_back(lh);
472 int geo = -1;
473 if (i < 8) {
474 geo = i * 2 + 1;
475 } else if (i < 20) {
476 geo = i * 2 + 2;
477 } else if (i < 36) {
478 geo = i * 2 + 3;
479 } else {
480 geo = i * 2 + 4;
481 }
482 lh->Gen(geo);
483
484 for(int j=0;j<WireNo[lyr]/2;j++){
485 /// wh is signal wire
486 MdcGeoWire *wh = new MdcGeoWire();
487 const int wireId = wircount+j;
488 wh->Id(wireId);
489
490 /// caluculate misaligned wire position
491 a1 = R[lyr]*cos(phi0+j*cellphi)+Sx[i_east]+wirePosVec[wireId][0];
492 a2 = R[lyr]*sin(phi0+j*cellphi)+Sy[i_east]+wirePosVec[wireId][1];
493 a3 = wLength[lyr]/2+wirePosVec[wireId][2]+Sz[i_east];
494 //a3 = wLength[lyr]/2+wirePosVec[wireId][2];
495 //MDC rotate around x firstly;
496 double xb = a1;
497 double yb = a2 * cos(Rx[i_east]) + a3 * sin(Rx[i_east]);
498 double zb = a3 * cos(Rx[i_east]) - a2 * sin(Rx[i_east]);
499
500 //then MDC rotate around y;
501 double xbnew = xb * cos(Ry[i_east]) - zb * sin(Ry[i_east]);
502 double ybnew = yb;
503 double zbnew = xb * sin(Ry[i_east]) + zb * cos(Ry[i_east]);
504 ///
505 noma1 = R[lyr]*cos(nomphi0+j*cellphi);
506 noma2 = R[lyr]*sin(nomphi0+j*cellphi);
507 /// attention wLength[lyr] not the true length of the wire
508 noma3 = wLength[lyr]/2;
509 //HepPoint3D pb(a1,a2,a3);
510 HepPoint3D pb(xbnew,ybnew,zbnew);
511 HepPoint3D nompb(noma1,noma2,noma3);
512 HepPoint3D wireposb(wirePosVec[wireId][0],wirePosVec[wireId][1],wirePosVec[wireId][2]);
513 wh->Backward(pb);
514 wh->nomBackward(nompb);
515 wh->BWirePos(wireposb);
516
517 //double g1 = noma1*cos(Rz[i_east]) - noma2*sin(Rz[i_east]) + Sx[i_east] + wirePosVec[wireId][0];
518 //double g2 = noma2*cos(Rz[i_east]) + noma1*sin(Rz[i_east]) + Sy[i_east] + wirePosVec[wireId][1];
519 //double g3 = noma3 + wirePosVec[wireId][2];
520 //HepPoint3D backward_test(g1,g2,g3);
521
522 a1 = R[lyr]*cos(phi0+(j+lh->Shift())*cellphi)+Sx[i_west]+wirePosVec[wireId][3];
523 a2 = R[lyr]*sin(phi0+(j+lh->Shift())*cellphi)+Sy[i_west]+wirePosVec[wireId][4];
524 //a3 = -wLength[lyr]/2+wirePosVec[wireId][5];
525 a3 = -wLength[lyr]/2+wirePosVec[wireId][5]+Sz[i_west];
526 //MDC rotate around x firstly;
527 double xf = a1;
528 double yf = a2 * cos(Rx[i_west]) + a3 * sin(Rx[i_west]);
529 double zf = a3 * cos(Rx[i_west]) - a2 * sin(Rx[i_west]);
530
531 //then MDC rotate around y;
532 double xfnew = xf * cos(Ry[i_west]) - zf * sin(Ry[i_west]);
533 double yfnew = yf;
534 double zfnew = xf * sin(Ry[i_west]) + zf * cos(Ry[i_west]);
535
536 noma1 = R[lyr]*cos(nomphi0+(j+lh->nomShift())*cellphi);
537 noma2 = R[lyr]*sin(nomphi0+(j+lh->nomShift())*cellphi);
538 noma3 = -wLength[lyr]/2;
539 //HepPoint3D pf(a1,a2,a3);
540 HepPoint3D pf(xfnew, yfnew, zfnew);
541 HepPoint3D nompf(noma1,noma2,noma3);
542 HepPoint3D wireposf(wirePosVec[wireId][3],wirePosVec[wireId][4],wirePosVec[wireId][5]);
543 wh->Forward(pf);
544 wh->nomForward(nompf);
545 wh->FWirePos(wireposf);
546
547 //double f1 = noma1*cos(Rz[i_west]) - noma2*sin(Rz[i_west]) + Sx[i_west] + wirePosVec[wireId][3];
548 //double f2 = noma2*cos(Rz[i_west]) + noma1*sin(Rz[i_west]) + Sy[i_west] + wirePosVec[wireId][4];
549 //double f3 = noma3 + wirePosVec[wireId][5];
550 //HepPoint3D forward_test(f1,f2,f3);
551 /// test
552 /// std::cout<<" wireId: "<<wireId<<" backward position:"<< wh->Backward()<<" forward position: "<<wh->Forward()<<std::endl;
553 wh->Layer(i);
554 wh->Cell(j);
555 wh->Stat(0);
556 wh->Slant(ebusing[lyr]);
557 wh->nomSlant(nomebusing[lyr]);
558 //cout<<ebusing[lyr]<<endl;
559 /// set the tension and wire sag
560 wh->Tension(wireTensionVec[wireId]);
561 //std::cout<<" wh->Tension: "<< wh->Tension()<<std::endl;
562 //std::cout<<" wire sag: "<<wh->Sag(wireId)<<std::endl;
563
564 /// because of just before have fLayers.push_back(lh);
565 /// so use it2 = fLayers.end()-1 to get the MdcGeoLayer
566 /// which have just push_back into the stack
567 vector<MdcGeoLayer*>::iterator it2 = fLayers.end()-1;
568 wh->Lyr(* it2);
569 fWires.push_back(wh);
570 }
571 }
572
573 fMisc.OuterR(fOutR[1]);
574 fMisc.InnerR(fInnerR[2]);
575 fMisc.OuterTk(fOutR[1]-fInnerR[1]);
576 fMisc.InnerTk(fOutR[2]-fInnerR[2]);
577 fMisc.NSWire(fWires.size());
578 fMisc.NFWire(TWireNo-fWires.size());
579
580 fMisc.LayerNo(TLayerNo);
581 fMisc.WireNo(TWireNo);
582 fMisc.SLayerNo(TSignalLayerNo);
583 fMisc.SWireR(signalWireR);
584 fMisc.FWireR(fieldWireR);
585
586 for(i=0;i<segmentNo;i++){
587 MdcGeoEnd * end=new MdcGeoEnd();
588 end->Id(i);
589 end->Length(fLength[i]);
590 end->InnerR(fInnerR[i]);
591 end->OutR(fOutR[i]);
592 end->Z(fZ[i]);
593 end->Name(fName[i]);
594 fEnd.push_back(end);
595 }
596}
597
598void MdcGeomSvc::ReadTensionDataBase(std::vector<double> & wireTensionVec){
599 std::string fullPath = "/Calib/MdcAlign";
600 MsgStream log(messageService(), name());
601 log << MSG::INFO<<"ReadTensionDataBase() wireTensionPath = "<<fullPath<< endreq;
602 cout << "Read wireTension from Calibration Database!" << endl;
603
604 SmartDataPtr<CalibData::MdcAlignData> tension(m_pCalibDataSvc, fullPath);
605 if( ! tension){
606 log << MSG::ERROR << "can not get MdcAlignConst via SmartPtr"
607 << endreq;
608 return;
609 }
610 for(int i=0;i<6796;i++){
611 double tens = tension->gettension(i);
612 wireTensionVec.push_back(tens);
613 }
614}
615void MdcGeomSvc::ReadWirePosDataBase(std::vector<vector<double> > & wirePosVec){
616 std::string fullPath = "/Calib/MdcAlign";
617 MsgStream log(messageService(), name());
618 log << MSG::INFO<<"ReadWirePosDataBase() wirePositionPath = "<<fullPath<< endreq;
619
620 cout << "Read wirePosition from Calibration Database!" << endl;
621 SmartDataPtr<CalibData::MdcAlignData> wirepos(m_pCalibDataSvc, fullPath);
622 if( ! wirepos){
623 log << MSG::ERROR << "can not get MdcAlignConst via SmartPtr"
624 << endreq;
625 return;
626 }
627 double dxe,dye,dze,dxw,dyw,dzw;
628 for(int j=0; j<6796; j++){
629
630 dxe = wirepos->getdxWireEast(j);
631 dye = wirepos->getdyWireEast(j);
632 dze = wirepos->getdzWireEast(j);
633 dxw = wirepos->getdxWireWest(j);
634 dyw = wirepos->getdyWireWest(j);
635 dzw = wirepos->getdzWireWest(j);
636
637 wirePosVec[j].push_back(dxe);
638 wirePosVec[j].push_back(dye);
639 wirePosVec[j].push_back(dze);
640 wirePosVec[j].push_back(dxw);
641 wirePosVec[j].push_back(dyw);
642 wirePosVec[j].push_back(dzw);
643 }
644
645}
646void MdcGeomSvc::ReadAliParDataBase(vector<double>& Sx, vector<double>& Sy, vector<double>& Sz,
647 vector<double>& Rx, vector<double>& Ry, vector<double>& Rz){
648 MsgStream log(messageService(), name());
649 std::string fullPath = "/Calib/MdcAlign";
650 log << MSG::INFO<<"ReadAliParDataBase() alignParPath = "<<fullPath<< endreq;
651 cout << "Read alignment parameters from Calibration Database!" << endl;
652
653 SmartDataPtr<CalibData::MdcAlignData> alignpar(m_pCalibDataSvc, fullPath);
654 if( ! alignpar){
655 log << MSG::ERROR << "can not get MdcAlignConst via SmartPtr"
656 << endreq;
657 return;
658 }
659 double tmp1,tmp2,tmp3,tmp4,tmp5,tmp6;
660 for(int k=0;k<16;k++){
661 tmp1 = alignpar->getdxEP(k);
662 tmp2 = alignpar->getdyEP(k);
663 tmp3 = alignpar->getdzEP(k);
664 tmp4 = alignpar->getrxEP(k);
665 tmp5 = alignpar->getryEP(k);
666 tmp6 = alignpar->getrzEP(k);
667
668 Sx.push_back(m_wholeShiftX+tmp1);
669 Sy.push_back(m_wholeShiftY+tmp2);
670 Sz.push_back(m_wholeShiftZ+tmp3);
671 Rx.push_back(m_wholeRotatX+tmp4);
672 Ry.push_back(m_wholeRotatY+tmp5);
673 Rz.push_back(m_wholeRotatZ+tmp6);
674 }
675
676}
677const int MdcGeomSvc::getWireSize()
678{
679 return fWires.size();
680}
681
682const int MdcGeomSvc::getLayerSize()
683{
684 return fLayers.size();
685}
686
687const int MdcGeomSvc::getSuperLayerSize()
688{
689 return fSupers.size();
690}
691
692const int MdcGeomSvc::getGeneralLayerSize()
693{
694 return fGenerals.size();
695}
696
697const int MdcGeomSvc::getSegmentNo()
698{
699 return fEnd.size();
700}
701
702void MdcGeomSvc::Dump(){}
703
704const int MdcGeomSvc::getAlignParIndexEast(int lyr) const{
705 /// the Mdc small endplate east
706 if((lyr>=0) && (lyr<=16)) return 0;
707 /// the 1th ring east
708 else if((lyr>=17) && (lyr<=20)) return 1;
709 /// the 2th ring east
710 else if((lyr>=21) && (lyr<=24)) return 2;
711 /// the 3th ring east
712 else if((lyr>=25) && (lyr<=28)) return 3;
713 /// the 4th ring east
714 else if((lyr>=29) && (lyr<=32)) return 4;
715 /// the 5th ring east
716 else if((lyr>=33) && (lyr<=36)) return 5;
717 /// the 6th ring east
718 else if((lyr>=37) && (lyr<=41)) return 6;
719 /// the Mdc big endplate east
720 else if(lyr>=42) return 7;
721 else std::cout<<" Hi, ERROR OCCUR !!!"<<std::endl;
722 return -1;
723}
724
725
726const int MdcGeomSvc::getAlignParIndexWest(int lyr) const{
727 /// the Mdc small endplate west
728 if((lyr>=0) && (lyr<=16)) return 8;
729 /// the 1th ring west
730 else if((lyr>=17) && (lyr<=20)) return 9;
731 /// the 2th ring west
732 else if((lyr>=21) && (lyr<=24)) return 10;
733 /// the 3th ring west
734 else if((lyr>=25) && (lyr<=28)) return 11;
735 /// the 4th ring west
736 else if((lyr>=29) && (lyr<=32)) return 12;
737 /// the 5th ring west
738 else if((lyr>=33) && (lyr<=36)) return 13;
739 /// the 6th ring west
740 else if((lyr>=37) && (lyr<=41)) return 14;
741 /// the Mdc big endplate west
742 else if(lyr>=42) return 15;
743 else std::cout<<" Hi, ERROR OCCUR !!!"<<std::endl;
744 return -1;
745}
746
747
748/// this handle function is prepared for special use
749void MdcGeomSvc::handle(const Incident& inc){
750 MsgStream log( messageService(), name() );
751 log << MSG::DEBUG << "handle: " << inc.type() << endreq;
752 IDataProviderSvc* m_eventSvc;
753 Gaudi::svcLocator()->service("EventDataSvc", m_eventSvc, true);
754 SmartDataPtr<Event::EventHeader> eventHeader(m_eventSvc,"/Event/EventHeader");
755 if (!eventHeader) {
756 log << MSG::FATAL << "Could not find Event Header" << endreq;
757 }
758 if (m_updataalign) return;
759 if (inc.type() == "NewRun" ){
760 log << MSG::DEBUG << "Begin Event" << endreq;
761 clean();
762 m_updataalign = true;
763 if(m_nomcalignment&&m_mindex==0) {
764 int RunNo=eventHeader->runNumber();
765 if(RunNo<0) m_readAlignParDataBase=false ;
766 else m_readAlignParDataBase=true;
767 m_mindex+=1;
768 cout<<"m__RunNo="<<RunNo<<"m_mindex="<<m_mindex<<endl;
769 }
770 //std::cout<<"############"<<m_readAlignParDataBase<<std::endl;
771 ReadFilePar();
772 AddWireNeighbors();
773 }
774}
775
776const MdcGeoWire * const
777 MdcGeomSvc::Wire(unsigned id){
778 if (id < fWires.size())
779 return fWires[id];
780
781 return 0;
782 }
783
784const MdcGeoWire * const
785 MdcGeomSvc::Wire(unsigned lyrid, unsigned wirid){
786 if ((lyrid <fLayers.size()) && ((int) wirid < Layer(lyrid)->NCell()))
787 return fWires[Layer(lyrid)->Wirst() + wirid];
788
789 return 0;
790 }
791
792const MdcGeoLayer * const
793 MdcGeomSvc::Layer(unsigned id){
794 if (id < fLayers.size())
795 return fLayers[id];
796
797 return 0;
798 }
799
800const MdcGeoSuper * const
801 MdcGeomSvc::SuperLayer(unsigned id){
802 if (id < fSupers.size())
803 return fSupers[id];
804
805 return 0;
806 }
807
808const MdcGeoGeneral * const
809 MdcGeomSvc::GeneralLayer(unsigned id){
810 if (id < fGenerals.size())
811 return & fGenerals[id];
812
813 return 0;
814 }
815
816const MdcGeoMisc * const
817MdcGeomSvc::Misc(void){
818 return & fMisc;
819}
820
821const MdcGeoEnd * const
822 MdcGeomSvc::End(unsigned id){
823 if (id < fEnd.size())
824 return fEnd[id];
825
826 return 0;
827 }
828
829
830bool MdcGeomSvc::getSagFlag(void) {
831
832 return m_doSag;
833}
834
835void MdcGeomSvc::AddWireNeighbors(){
836 std::vector<double> BackwardPhi;
837 std::vector<double> ForwardPhi;
838 std::vector<double> Center_r;
839 std::vector<double> Cho;//half distance in r-phi plane betweem two end points;
840 int layerID = -1;
841 MdcGeoWire* WireK ;
842 for(unsigned k = 0;k<fWires.size();k++) {
843 WireK = fWires.at(k);
844 if(layerID!=WireK->Layer()){
845 Center_r.push_back(sqrt(pow(WireK->Backward().x()/2.+WireK->Forward().x()/2.,2.)+pow(WireK->Backward().y()/2.+WireK->Forward().y()/2.,2.)));
846 Cho.push_back(0.5*sqrt(pow(WireK->Backward().x()-WireK->Forward().x(),2.)+pow(WireK->Backward().y()-WireK->Forward().y(),2.)));
847 layerID=WireK->Layer();
848 }
849 double phi1=0;
850 double phi2=0;
851 if(WireK->Backward().y()>=0){
852 phi1 = acos(WireK->Backward().x()/sqrt(pow(WireK->Backward().x(),2.)+pow(WireK->Backward().y(),2.)));
853 }
854 else {
855 phi1 =2*M_PI-acos(WireK->Backward().x()/sqrt(pow(WireK->Backward().x(),2.)+pow(WireK->Backward().y(),2.)));
856 }
857
858 if(WireK->Forward().y()>=0){
859 phi2 = acos(WireK->Forward().x()/sqrt(pow(WireK->Forward().x(),2.)+pow(WireK->Forward().y(),2.)));
860 }
861 else {
862 phi2 =2*M_PI-acos(WireK->Forward().x()/sqrt(pow(WireK->Forward().x(),2.)+pow(WireK->Forward().y(),2.)));
863 }
864
865 BackwardPhi.push_back(phi1);
866 ForwardPhi.push_back(phi2);
867
868 }
869 for(unsigned k = 0;k<fWires.size();k++) {
870 int wireID;
871 int neiborType;
872 double dPhi,nPhi;
873 double phiBmax,phiBmin,phiFmax,phiFmin;
874 int layer1;
875 double Bphi,Fphi,BphiN,FphiN;
876 MdcGeoWire* WireN;
877 Bphi = BackwardPhi.at(k);
878 Fphi = ForwardPhi.at(k);
879 WireK = fWires.at(k);
880 if(k==0||fWires.at(k-1)->Layer()!=WireK->Layer()){
881 wireID = k+WireK->Lyr()->NCell()-1;
882 neiborType =1;
883 WireK->MdcGeoWire::AddWireNeighbor(fWires.at(wireID),wireID,neiborType);
884 wireID = k+1;
885 neiborType =2;
886 WireK->MdcGeoWire::AddWireNeighbor(fWires.at(wireID),wireID,neiborType);
887 }
888 else if(k==6795||fWires.at(k+1)->Layer()!=WireK->Layer()){
889 wireID = k-1;
890 neiborType =1;
891 WireK->MdcGeoWire::AddWireNeighbor(fWires.at(wireID),wireID,neiborType);
892 wireID = k-WireK->Lyr()->NCell()+1;
893 neiborType =2;
894 WireK->MdcGeoWire::AddWireNeighbor(fWires.at(wireID),wireID,neiborType);
895 }
896 else{
897 wireID = k-1;
898 neiborType =1;
899 WireK->MdcGeoWire::AddWireNeighbor(fWires.at(wireID),wireID,neiborType);
900 wireID = k+1;
901 neiborType =2;
902 WireK->MdcGeoWire::AddWireNeighbor(fWires.at(wireID),wireID,neiborType);
903 }
904 for(int n=k-1;n>=0&&(WireK->Layer()-fWires.at(n)->Layer())<2;n--){
905 int flag=0;
906 WireN = fWires.at(n);
907 BphiN = BackwardPhi.at(n);
908 FphiN = ForwardPhi.at(n);
909 if((WireK->Layer()-WireN->Layer())==1){
910 layer1 = WireK->Layer();
911 dPhi = (M_PI/WireN->Lyr()->NCell() + M_PI/WireK->Lyr()->NCell())*m_NeighborRange ;
912 nPhi = atan(Cho.at(layer1)/Center_r.at(layer1))-atan((WireN->Forward().z()/WireK->Forward().z())*Cho.at(layer1)/Center_r.at(layer1));//angular offset beacause of different wire lenghth
913 if((Bphi>Fphi&&(Bphi-Fphi)<M_PI)||(Fphi-Bphi)>M_PI){
914 phiBmax = CalculatePhi(Bphi+dPhi-nPhi);
915 phiBmin = CalculatePhi(Bphi-dPhi-nPhi);
916 phiFmax = CalculatePhi(Fphi+dPhi+nPhi);
917 phiFmin = CalculatePhi(Fphi-dPhi+nPhi);
918 }
919 else{
920 phiBmax = CalculatePhi(Bphi+dPhi+nPhi);
921 phiBmin = CalculatePhi(Bphi-dPhi+nPhi);
922 phiFmax = CalculatePhi(Fphi+dPhi-nPhi);
923 phiFmin = CalculatePhi(Fphi-dPhi-nPhi);
924 }
925 wireID = n;
926 neiborType =3;
927 if(((BphiN>phiBmin&&BphiN<phiBmin+M_PI/2.)||BphiN<phiBmin-3*M_PI/2.)&&((FphiN<phiFmax&&FphiN>phiFmax-M_PI/2.)||FphiN>phiFmax+3*M_PI/2.)){
928 WireK->MdcGeoWire::AddWireNeighbor(WireN,wireID,neiborType);
929 flag=1;
930 }
931 if(((FphiN>phiFmin&&FphiN<phiFmin+M_PI/2.)||FphiN<phiFmin-3*M_PI/2.)&&((BphiN<phiBmax&&BphiN>phiBmax-M_PI/2.)||BphiN>phiBmax+3*M_PI/2.)){
932 if(flag==0) WireK->MdcGeoWire::AddWireNeighbor(WireN,wireID,neiborType);
933 }
934 }
935 }
936 for(int n=k+1;n<6796&&(fWires.at(n)->Layer()-WireK->Layer())<2;n++){
937 int flag=0;
938 WireN = fWires.at(n);
939 BphiN = BackwardPhi.at(n);
940 FphiN = ForwardPhi.at(n);
941 if((WireN->Layer()-WireK->Layer())==1){
942 layer1 = WireN->Layer();
943 dPhi = (M_PI/WireK->Lyr()->NCell() + M_PI/WireN->Lyr()->NCell())*m_NeighborRange;
944 nPhi = atan(Cho.at(layer1)/Center_r.at(layer1))-atan((WireK->Forward().z()/WireN->Forward().z())*Cho.at(layer1)/Center_r.at(layer1));
945 if((BphiN>FphiN&&(BphiN-FphiN)<M_PI)||(FphiN-BphiN)>M_PI){
946 phiBmax = CalculatePhi(BphiN+dPhi-nPhi);
947 phiBmin = CalculatePhi(BphiN-dPhi-nPhi);
948 phiFmax = CalculatePhi(FphiN+dPhi+nPhi);
949 phiFmin = CalculatePhi(FphiN-dPhi+nPhi);
950 }
951 else{
952 phiBmax = CalculatePhi(BphiN+dPhi+nPhi);
953 phiBmin = CalculatePhi(BphiN-dPhi+nPhi);
954 phiFmax = CalculatePhi(FphiN+dPhi-nPhi);
955 phiFmin = CalculatePhi(FphiN-dPhi-nPhi);
956 }
957 wireID = n;
958 neiborType =4;
959 if(((Bphi>phiBmin&&Bphi<phiBmin+M_PI/2.)||Bphi<phiBmin-3*M_PI/2.)&&((Fphi<phiFmax&&Fphi>phiFmax-M_PI/2.)||Fphi>phiFmax+3*M_PI/2.)){
960 WireK->MdcGeoWire::AddWireNeighbor(WireN,wireID,neiborType);
961 flag=1;
962 }
963 if(((Fphi>phiFmin&&Fphi<phiFmin+M_PI/2.)||Fphi<phiFmin-3*M_PI/2.)&&((Bphi<phiBmax&&Bphi>phiBmax-M_PI/2.)||Bphi>phiBmax+3*M_PI/2.)){
964 if(flag==0) WireK->MdcGeoWire::AddWireNeighbor(WireN,wireID,neiborType);
965 }
966 }
967 }
968
969 }
970}
971
972void MdcGeomSvc::RenewWireNeighbors(double range){
973 m_NeighborRange = range;
974 for(unsigned k = 0;k<fWires.size();k++) {
975 fWires.at(k)->ClearNeighbors();
976 }
977 MdcGeomSvc::AddWireNeighbors();
978}
979
980 const double MdcGeomSvc::CalculatePhi(double phi) const{
981 if(0<phi&&phi<2*M_PI){
982 return phi;
983 }
984 else if(phi>2*M_PI){
985 return phi-2*M_PI;
986 }
987 else {
988 return phi+2*M_PI;
989 }
990 }
991
992
993
994/**
995 * @brief Get the Intersection Point, assuming line A, B are in same surface,
996 * and A1 , B1 represents a point on line A,B.
997 * dA, dB is direction vector of A and B. that is, dA can be A2-A1 where A1, A2 is 2 points on line A
998 * line A, B is not parallel (they meets).
999 * @return HepPoint3D
1000 */
1001static HepPoint3D getIntersection(HepPoint3D A1,HepPoint3D dA, HepPoint3D B1, HepPoint3D dB){
1002 // making them in direction z+ and unit in z
1003 dA /= dA.z();
1004 dB /= dB.z();
1005 // moving b1 same z as A in the direction
1006 B1 += dB *(A1.z()-B1.z());
1007
1008 HepPoint3D deltaPoint=B1-A1;
1009 HepPoint3D deltaDir=dB-dA;
1010 double deltaZ=-deltaDir.dot(deltaPoint)/deltaDir.mag2();
1011 return A1 + deltaZ* dA;
1012
1013}
1014
1015static HepPoint3D interpolate(HepPoint3D A1,HepPoint3D A2, double targetZ){
1016 HepPoint3D dA=A2 - A1;
1017
1018 return A1 + dA*( (targetZ - A1.z()) / dA.z() );
1019
1020}
1021static double interpolate(double A1, double Za1,double A2,double Za2, double targetZ){
1022 double dA=A2 - A1;
1023 double dz=Za2-Za1;
1024
1025 return A1 + dA*( (targetZ - Za1) / dz );
1026
1027}
1028
1029static double squared(double x){return x*x;}
1030
1031
1032
1033/**
1034 * @brief get touching shape of 2 adjacent layer wires.
1035 * return vertexs of intersection shape if the wires' sensitive area connects.
1036 * return empty if not.
1037 * points does go in order
1038
1039 *
1040 * @param wireA
1041 * @param wireB
1042 * @return vector<HepPoint3D> return 4 vertexs of intersection parallelogram if the wires' sensitive area connects.
1043 * return empty if not.
1044 */
1045vector<HepPoint3D> MdcGeomSvc::getAdjacentIntersectionShape(const MdcGeoWire* wireA,const MdcGeoWire* wireB) {
1046 vector<HepPoint3D> ret;
1047
1048 //to dbg: b MdcGeomSvc.cxx:MdcGeomSvc::getAdjacentIntersectionShape
1049
1050 HepPoint3D A1=wireA->Forward();
1051 HepPoint3D A2=wireA->Backward(); // OOps!!!
1052 HepPoint3D B1=wireB->Forward(); // is -
1053 HepPoint3D B2=wireB->Backward(); // is +
1054
1055 double maxZAllowed=(
1056 std::max(wireA->Backward().z(),wireB->Backward().z())
1057 // std::max(wireA->Forward().z(), wireB->Forward().z())
1058 );
1059 double minZAllowed=(
1060 std::min(wireA->Forward().z(), wireB->Forward().z())
1061 // std::min(wireA->Backward().z(),wireB->Backward().z())
1062 );
1063
1064 HepPoint3D directionR=(A1+A2+B1+B2)/4;// in general, should be (almost) perpendicular to line A and B
1065 directionR.setZ(0);
1066 directionR.setMag(1); // now a directional vector.
1067 // HepPoint3D directionZ(0,0,1);
1068 // HepPoint3D directionPhi=directionR.cross(directionZ);
1069
1070 const MdcGeoWire* wireAl=this->Wire(static_cast<unsigned int>(wireA->GetNeighborIDType1()));
1071 //const MdcGeoWire* wireA2=this->Wire(static_cast<unsigned int>(wireA->GetNeighborIDType2()));
1072 const MdcGeoWire* wireB1=this->Wire(static_cast<unsigned int>(wireB->GetNeighborIDType1()));
1073 // const MdcGeoWire* wireB2=this->Wire(static_cast<unsigned int>(wireB->GetNeighborIDType2()));
1074 double typicalWidthA=(wireAl->Forward() - wireA->Forward()).mag();
1075 double typicalWidthB=(wireB1->Forward() - wireB->Forward()).mag(); // using dist at endplate;
1076
1077 // double radiusLayerA=this->Layer(wireA->Layer())->Radius();
1078 // double radiusLayerB=this->Layer(wireB->Layer())->Radius();
1079 // double typicalLayerThickness=(radiusLayerB-radiusLayerA)/(wireB->Layer()-wireA->Layer());
1080 // double dWAdr=typicalWidthA/typicalLayerThickness*2; // hope we have a better alternative.
1081 // double dWBdr=typicalWidthB/typicalLayerThickness*2;
1082
1083 // if we have r implemented, this will be calced from dr and dWAdr
1084 double widthAtSurfaceA=typicalWidthA;
1085 double widthAtSurfaceB=typicalWidthB;
1086 // double minWidth=min(widthAtSurfaceA,widthAtSurfaceB);
1087
1088 //double centerZ=getClosestZ(wireA,wireB);
1089 //double rA=wireA->GetX(centerZ);
1090
1091 double surfaceA=directionR.dot((A1+A2)/2); // centerZ
1092 double surfaceB=directionR.dot((B1+B2)/2);
1093 // double surfaceMid=(surfaceA+surfaceB)/2;
1094 double surfaceDiff=(surfaceB-surfaceA);
1095
1096 // so the midsurface: any X that directionR.dot(X)==surfaceMid
1097 HepPoint3D A1m=A1+ surfaceDiff/2*directionR;
1098 HepPoint3D A2m=A2+ surfaceDiff/2*directionR;
1099 HepPoint3D B1m=B1- surfaceDiff/2*directionR;
1100 HepPoint3D B2m=B2- surfaceDiff/2*directionR;
1101 // HepPoint3D d1m=
1102
1103 HepPoint3D direction1m=B1m - A1m;
1104 direction1m.setZ(0);
1105 {
1106 HepPoint3D direction2m=B2m - A2m;
1107 direction2m.setZ(0);
1108 if (direction2m.mag2() > direction1m.mag()) direction1m = direction2m;
1109 }
1110 direction1m.setMag(1);
1111
1112 //HepPoint3D center=getIntersection(A1m,A2m, B1m, B2m );
1113 HepPoint3D directionA=(A2-A1);directionA/=directionA.z();
1114 HepPoint3D directionB=(B2-B1);directionB/=directionB.z();
1115 HepPoint3D deltaA=direction1m*(widthAtSurfaceA/2); // A1->B1 direction
1116 HepPoint3D deltaB=direction1m*(widthAtSurfaceB/2); // A1->B1 direction
1117
1118 if ((directionA-directionB).mag2()<= 0.01*0.01){ // we know wire distance ~ 10 mm so a 10^-3 level error is allowed.
1119 // if parallel
1120 if (direction1m.mag2()>squared(widthAtSurfaceA+widthAtSurfaceB)/4) {// have no relation at all.
1121 return ret;
1122 // donothing. let it return.
1123 }
1124 else{// has interection
1125 ret.push_back(B1m-deltaB);
1126 ret.push_back(A1m+deltaA);
1127 ret.push_back(A2m+deltaA);
1128 ret.push_back(B2m-deltaB);
1129 }// if parallel and has intersection: intersection is a parallelogram, with a pair of lines tangent to Z
1130 }
1131 else{// expect a parallelogram. Z: not now. we want 3, 4, 5 or 6 angle shape . take care!
1132
1133 HepPoint3D Ptop=getIntersection(A1m+deltaA,directionA,B1m-deltaB,directionB); // appears wrong size..
1134 HepPoint3D Prht=getIntersection(A1m+deltaA,directionA,B1m+deltaB,directionB);
1135 HepPoint3D Plft=getIntersection(A1m-deltaA,directionA,B1m-deltaB,directionB);
1136 HepPoint3D Pbot=getIntersection(A1m-deltaA,directionA,B1m+deltaB,directionB);
1137
1138 if (Ptop.z() < Pbot.z()) std::swap(Ptop,Pbot);
1139
1140 HepPoint3D Pmtp=Prht.z()>Plft.z()? Prht: Plft; // median top
1141 HepPoint3D Pmbt=Prht.z()>Plft.z()? Plft: Prht; // median bottom
1142
1143 HepPoint3D PointsAtMaxZ[2]; //maxz or top
1144 int PointsAtMaxZCount=0;
1145 HepPoint3D PointsAtMinZ[2];
1146 int PointsAtMinZCount=0;
1147 HepPoint3D PointsLeftGood[1]={Plft};
1148 int PointsLeftGoodCount=1;
1149 HepPoint3D PointsRightGood[1]={Prht};
1150 int PointsRightGoodCount=1;
1151
1152 // a chain of ifs
1153 if (Ptop.z()<=maxZAllowed){
1154 PointsAtMaxZ[0]=Ptop;
1155 PointsAtMaxZCount=1;
1156 }else if (Ptop.z()>maxZAllowed && maxZAllowed>=Pmtp.z() ){ // only this 2
1157 PointsAtMaxZ[0]=interpolate(Ptop,Plft,maxZAllowed);//L
1158 PointsAtMaxZ[1]=interpolate(Ptop,Prht,maxZAllowed);//R
1159 PointsAtMaxZCount=2;
1160 }else if (Pmtp.z()>maxZAllowed && maxZAllowed>= Pmbt.z() ){ // midian
1161 PointsAtMaxZ[0]=interpolate(Ptop,Pmbt,maxZAllowed);// Top-mbot
1162 PointsAtMaxZ[1]=interpolate(Pmtp,Pbot,maxZAllowed);// mrop-Bot
1163 PointsAtMaxZCount=2;
1164 if (Plft==Pmtp) {// lft is higher
1165 std::swap(PointsAtMaxZ[0],PointsAtMaxZ[1]); // now 0=Left, 0-Right
1166 PointsLeftGoodCount=0; // so remove left point
1167 }else{ // rht is higher
1168 PointsRightGoodCount=0;
1169 }
1170
1171 // glad i choose not to play it in order
1172 }else if (Pmbt.z()>maxZAllowed && maxZAllowed>= Pbot.z() ){ // lower
1173 PointsAtMaxZ[0]=interpolate(Pbot,Plft,maxZAllowed);//L
1174 PointsAtMaxZ[1]=interpolate(Pbot,Prht,maxZAllowed);//R
1175 PointsAtMaxZCount=2;
1176 PointsLeftGoodCount=PointsRightGoodCount=0; // remove both
1177 }else{ // no---- we got nothing at all
1178 return ret;
1179 }
1180 // and then with minZAllowed
1181 if ( minZAllowed <=Pbot.z() ){ //
1182 PointsAtMinZ[0]=Pbot;
1183 PointsAtMinZCount=1;
1184 }else if (Pbot.z() < minZAllowed && minZAllowed <=Pmbt.z()){ // lower.
1185 PointsAtMinZ[0]=interpolate(Pbot,Plft,minZAllowed);//L //TODO: change to Pbot
1186 PointsAtMinZ[1]=interpolate(Pbot,Prht,minZAllowed);//R //TODO: change to Pbot
1187 PointsAtMinZCount=2;
1188 }else if (Pmbt.z() < minZAllowed && minZAllowed <=Pmtp.z()){// mid.
1189 PointsAtMinZ[0]=interpolate(Ptop,Pmbt,minZAllowed);// Top-mbot
1190 PointsAtMinZ[1]=interpolate(Pmtp,Pbot,minZAllowed);// mrop-Bot
1191 PointsAtMinZCount=2;
1192 if (Plft==Pmtp) {// lft is higher
1193 std::swap(PointsAtMaxZ[0],PointsAtMaxZ[1]); // now 0=Left, 0-Right
1194 PointsRightGoodCount=0; // so remove lower point
1195 }else{ // rht is higher
1196 PointsLeftGoodCount=0;
1197 }
1198 }else if (Pmtp.z() < minZAllowed && minZAllowed <=Ptop.z()){// higher
1199 PointsAtMinZ[0]=interpolate(Ptop,Plft,minZAllowed);//L
1200 PointsAtMinZ[1]=interpolate(Ptop,Prht,minZAllowed);//R
1201 PointsAtMinZCount=2;
1202 PointsRightGoodCount=PointsLeftGoodCount=0;
1203 }else{
1204 return ret;
1205 }
1206
1207 // now insert items in clockwise order.
1208 for (int i=0;i<PointsAtMaxZCount;i++){
1209 ret.push_back(PointsAtMaxZ[i]);
1210 }
1211 for (int i=0;i<PointsRightGoodCount;i++){
1212 ret.push_back(PointsRightGood[i]);
1213 }
1214 for (int i=PointsAtMinZCount-1;i>=0;i--){//
1215 ret.push_back(PointsAtMinZ[i]);
1216 }
1217 for (int i=0;i<PointsLeftGoodCount;i++){
1218 ret.push_back(PointsLeftGood[i]);
1219 }
1220
1221
1222 }
1223 return ret;
1224
1225}
1226/**
1227 * @brief get 2 points on circle, phi1 and phi2, dist<pi by moving one of point -2pi to allow correct cross-zero interpolation, etc.
1228 *
1229 * return bool: had i moved a point?
1230 */
1231static bool circularPhiExtend(double& phi1, double& phi2){
1232 double& phiU= (phi1>phi2 ) ?phi1:phi2; // upper
1233 double& phiL= (phi1>phi2 ) ?phi2:phi1; // lower
1234 // a more (appeared) general method: make the phi dist < pi
1235 if (phiU-phiL<M_PI){
1236 phiU-=M_PI*2;
1237 return true;
1238 }
1239 return false;
1240}
1241
1242/**
1243 * @brief find the z where points on wireA, wireB has same phi. point is given as avg(A(z),B(z))
1244 *
1245 * @param wireA
1246 * @param wireB
1247 * @return HepPoint3D
1248 */
1249HepPoint3D MdcGeomSvc::getAdjacentIntersectionPointFast(const MdcGeoWire* wireA,const MdcGeoWire* wireB){
1250 HepPoint3D A1=wireA->Forward();
1251 HepPoint3D A2=wireA->Backward(); // OOps!!!
1252 HepPoint3D B1=wireB->Forward(); // is -
1253 HepPoint3D B2=wireB->Backward(); // is +
1254 double phiA1=A1.getPhi();
1255 double phiA2=A2.getPhi();
1256 double phiB1=B1.getPhi();
1257 double phiB2=B2.getPhi();
1258 circularPhiExtend(phiA1, phiA2);
1259 circularPhiExtend(phiB1, phiB2);
1260 double dPhiA=phiA2-phiA1;
1261 double dZA=A2.z()-A1.z();
1262 double phiAatB1=phiA1+(B1.z()-A1.z())/dZA*dPhiA;
1263 double dPhiB=phiB2-phiB1;
1264 double dZB=B2.z()-B1.z();
1265 // we want z that makes phiA(z)==phiB(z),
1266 // that is phiAatB1+(z-zB1)*dPhiA/dZA == phiB1+(z-zB1)*dPhiB/dZB
1267 double z=B1.z()+(phiAatB1-phiB1)/(dPhiB/dZB-dPhiA/dZA);
1268 return HepPoint3D(
1269 (wireA->GetX(z)+wireB->GetX(z))/2,
1270 (wireA->GetY(z)+wireB->GetY(z))/2,
1271 z
1272 );
1273}
1274
1275HepPoint3D MdcGeomSvc::getAdjacentIntersectionPointSlower(const MdcGeoWire* wireA,const MdcGeoWire* wireB){
1276 vector<double> bounds = getAdjacentIntersectionBounds(wireA,wireB);
1277 HepPoint3D ret(1./0,1./0,1./0);
1278 if (bounds.empty()) return ret;
1279 double z=(bounds[2]+bounds[3])/2;
1280 return HepPoint3D(
1281 (wireA->GetX(z)+wireB->GetX(z))/2,
1282 (wireA->GetY(z)+wireB->GetY(z))/2,
1283 z
1284 );
1285
1286}
1287
1288// circular problem:
1289
1290/**
1291 * @brief get bounding box in \varphi and z of touching shape of 2 adjacent layer wires.
1292 *
1293 * @param wireA
1294 * @param wireB
1295 * @return vector<double> empty if no touch. [phiL,phiU, zL,zU] if touching.
1296 * may have phiL<0 if the box goes over phi=0. (circular)
1297 */
1298vector<double> MdcGeomSvc::getAdjacentIntersectionBounds(const MdcGeoWire* wireA,const MdcGeoWire* wireB){
1299 vector<double> ret;
1300 vector<HepPoint3D> vertexs=getAdjacentIntersectionShape(wireA, wireB);
1301 if (vertexs.size()==0) return ret;
1302 else{
1303 // double maxZAllowed=std::max(
1304 // std::max(wireA->Backward().z(),wireB->Backward().z()),
1305 // std::max(wireA->Forward().z(), wireB->Forward().z())
1306 // );
1307 // double minZAllowed=std::min(
1308 // std::min(wireA->Forward().z(), wireB->Forward().z()),
1309 // std::min(wireA->Backward().z(),wireB->Backward().z())
1310 // );
1311
1312 // vertexs.push_back(vertexs[0]);
1313
1314 double phiL,phiU, zL,zU;
1315 bool firstp=true;
1316
1317 for (size_t i=0;i<vertexs.size()-1;i++){
1318 HepPoint3D p1=vertexs[i];
1319 // HepPoint3D p2=vertexs[i+1];
1320 // int sgn1=0;
1321 // int sgn2=0; // for vertex 1,2. 0 for in range, 1 for >max, -1 for <min
1322 // if (p1.z() > maxZAllowed) sgn1=1;
1323 // else if (p1.z() < minZAllowed) sgn1 = -1;
1324 // if (p2.z() > maxZAllowed) sgn2=1;
1325 // else if (p2.z() < minZAllowed) sgn2 = -1;
1326
1327 // if (sgn1==sgn2 && sgn1 !=0){// both at outside
1328 // continue;
1329 // }
1330 // else if (sgn1==sgn2 && sgn1 ==0){ // both at inside
1331
1332 // }else {// crossed angle
1333 // HepPoint3D pu,pl;
1334 // pu = interpolate(p1,p2,maxZAllowed);
1335 // pl = interpolate(p1,p2,minZAllowed);
1336 // if (sgn1 ==1) p1 = pu;
1337 // if (sgn1 ==-1) p1 = pl;
1338 // if (sgn2 ==1) p2 = pu;
1339 // if (sgn2 ==-1) p2 = pl;
1340
1341 // }
1342
1343 if (firstp){ // hate this copy-pasting
1344 HepPoint3D vertex=p1;
1345 zL=vertex.z();
1346 zU=vertex.z();
1347 phiL=vertex.getPhi();
1348 phiU=vertex.getPhi();
1349 firstp=false;
1350 }
1351 else{
1352 HepPoint3D vertex=p1;
1353 zL=std::min(zL,vertex.z());
1354 zU=std::max(zU,vertex.z());
1355 phiL=std::min(phiL,vertex.getPhi());
1356 phiU=std::max(phiU,vertex.getPhi());
1357 }
1358 }
1359
1360
1361 if (phiU >=M_PI *1.5 && phiL <= M_PI *.5) {//damn hated circular phi problem
1362 HepPoint3D vertex=vertexs[0];
1363 double phi2=vertex.getPhi();
1364 if (phi2>M_PI *1.5)phi2 -=M_PI *2;
1365 phiL=phiU=phi2;
1366 for (size_t i=1;i<vertexs.size();i++){
1367 HepPoint3D vertex=vertexs[i];
1368 double phi2=vertex.getPhi();
1369 if (phi2>M_PI *1.5)phi2 -=M_PI *2;
1370 phiL=std::min(phiL,phi2);
1371 phiU=std::max(phiU,phi2);
1372 }
1373 }
1374
1375 ret.push_back(phiL);
1376 ret.push_back(phiU);
1377 ret.push_back(zL);
1378 ret.push_back(zU);
1379 return ret;
1380 }
1381}
1382
1383
1384void MdcGeomSvc::testAdjacentIntersection(){
1385 clock_t time0=clock();
1386 ofstream fileout("testAdjacentIntersection.csv");
1387 MsgStream log(messageService(), name());
1388 log << MSG::WARNING<<__PRETTY_FUNCTION__<<endreq;
1389 const char * comma=", ";
1390 const char * endll="\n";
1391 fileout << "# A_layer" << comma
1392 << "A_wireID" << comma
1393 << "A_WirePhi(0)" << comma
1394 // << "A_FWirePos_phi" << comma
1395 // << "A_BWirePos_phi" << comma
1396 << "A_Forward_phi" << comma
1397 << "A_Backward_phi" << comma
1398 << "B_layer" << comma
1399 << "B_wireID" << comma
1400 << "B_WirePhi(0)" << comma
1401 << "B_Forward_phi" << comma
1402 << "B_Backward_phi" << comma
1403
1404 << "phiL,phiU, zL,zU" << comma
1405
1406 << "wireA.x0, wireA.y0"<<comma
1407 << "wireB.x0, wireB.y0"<<comma
1408
1409 << "line.a,line.b,line.c,line.x0, line.y0, line.z0, line.Angle2d, line.angle3d, line.cosAngleR"<< comma
1410 << "line.a,line.b,line.c,line.x0, line.y0, line.z0, line.Angle2d, line.angle3d, line.cosAngleR"<< comma
1411 << "line.a,line.b,line.c,line.x0, line.y0, line.z0, line.Angle2d, line.angle3d, line.cosAngleR"<< comma
1412 << "line.a,line.b,line.c,line.x0, line.y0, line.z0, line.Angle2d, line.angle3d, line.cosAngleR"<< comma
1413
1414 << "xa1,ya1,za1" << comma // wireA.Forward
1415 << "xa2,ya2,za2" << comma // wireA.Backward
1416 << "xb1,yb1,zb1" << comma // wireB.Forward
1417 << "xb2,yb2,zb2" << comma // wireB.Backward
1418 << "x1,y1,z1" << comma //getAdjacentIntersectionShape.x,y,z
1419 << "x2,y2,z2" << comma
1420 << "x3,y3,z3" << comma
1421 << "x4,y4,z4" << endll;
1422
1423 fileout << "# xn yn zn are coordinates of vertex" << endll;
1424 fileout << "# layer: ";
1425 for (unsigned i=0;i<43;++i){
1426 fileout<<comma<<i;
1427 }fileout << endll;
1428 fileout << "# layerRadius: ";
1429 for (unsigned i=0;i<43;++i){
1430 double r=this->Layer(i)->Radius();
1431 fileout<<comma<<r;
1432 }fileout << endll;
1433 fileout << "# distToNextLayerRadius: ";
1434 for (unsigned i=0;i<42;++i){
1435 double r=this->Layer(i+1)->Radius()-this->Layer(i)->Radius();
1436 fileout<<comma<<r;
1437 }fileout << endll;
1438
1439 for (size_t _wire_index=0;_wire_index<fWires.size();++_wire_index){
1440 const MdcGeoWire * wire=fWires[_wire_index];
1441 // if (wire->Id()!=1726) continue;
1442 std::vector<MdcGeoWire *> cross_neighbors=wire->GetNeighborType4();
1443 for (size_t _neighbor_index=0;_neighbor_index<cross_neighbors.size();_neighbor_index++){
1444 const MdcGeoWire * neighbor=cross_neighbors[_neighbor_index];
1445 if (abs(neighbor->Layer()-wire->Layer())!=1) continue;
1446 std::vector<HepPoint3D> shape=getAdjacentIntersectionShape(wire, neighbor);
1447 std::vector<double> bound=getAdjacentIntersectionBounds(wire, neighbor);
1448 vector<AdjCandiTgtLine2D> psblines=getAdjacentCandidateWithDriftCircle(wire,neighbor,10,10);
1449 fileout << wire->Layer() << comma
1450 << wire->Id() << comma
1451 << wire->WirePhi(0) << comma
1452 // << wire->FWirePos().getPhi()<< comma
1453 // << wire->BWirePos().getPhi()<< comma
1454 << wire->Forward().getPhi() << comma
1455 << wire->Backward().getPhi() << comma
1456
1457 << neighbor->Layer() << comma
1458 << neighbor->Id() << comma
1459 << neighbor->WirePhi(0) << comma
1460 << neighbor->Forward().getPhi() << comma
1461 << neighbor->Backward().getPhi() << comma;
1462 // << neighbor->FWirePos().getPhi()<<comma
1463 // << neighbor->BWirePos().getPhi()<<comma;
1464 // fileout<<"seg"<<comma;
1465 if (bound.empty()) {
1466 for (size_t i = 0; i < 4; ++i) fileout << 0. << comma;
1467 } else {
1468 for (size_t i = 0; i < bound.size(); i++){
1469 fileout << bound.at(i) << comma;
1470 }
1471 }
1472 // fileout<<"seg"<<comma;
1473 //"wireA.x0, wireA.y0, wireB.x0, wireB.y0"
1474 double Z0=0;
1475 if (!psblines.empty())Z0=psblines[0].z0;
1476
1477 fileout
1478 << wire ->GetX(Z0)<<comma
1479 << wire ->GetY(Z0)<<comma
1480 << neighbor->GetX(Z0) << comma
1481 << neighbor->GetY(Z0) << comma;
1482 // fileout<<"seg"<<comma;
1483 {
1484 for (size_t i = 0; i < psblines.size(); ++i)
1485 fileout
1486 <<psblines[i].a << comma
1487 <<psblines[i].b << comma
1488 <<psblines[i].c << comma
1489 <<psblines[i].x0 << comma
1490 <<psblines[i].y0 << comma
1491 <<psblines[i].z0 << comma
1492 <<psblines[i].angle2D << comma
1493 <<psblines[i].angle3D << comma
1494 <<psblines[i].angleWithR() << comma;
1495 for (size_t i = psblines.size(); i < 4; ++i)
1496 for (size_t ii = 0; ii < 9; ii++)
1497 fileout<<'0'<<comma;
1498 }
1499 // fileout<<"seg"<<comma;
1500 {
1501 HepPoint3D vertex;
1502 vertex=wire->Forward();
1503 fileout << vertex.x() << comma << vertex.y() << comma << vertex.z() <<comma;
1504 vertex=wire->Backward();
1505 fileout << vertex.x() << comma << vertex.y() << comma << vertex.z() <<comma;
1506 vertex=neighbor->Forward();
1507 fileout << vertex.x() << comma << vertex.y() << comma << vertex.z() <<comma;
1508 vertex=neighbor->Backward();
1509 fileout << vertex.x() << comma << vertex.y() << comma << vertex.z();// <<comma;
1510 }
1511 // fileout<<comma<<"seg";
1512 // if (shape.empty) // i dont want it to output. a csv don't need to be full-length
1513 for (size_t i = 0; i < shape.size(); i++){
1514 HepPoint3D vertex=shape[i];
1515 fileout <<comma<< vertex.x() << comma << vertex.y() << comma << vertex.z();
1516 //if (i+1<shape.size()) fileout<<comma;
1517 // alittle fix when shape is empty and it will lead to a trailing comma.
1518 }
1519 fileout<< endll;
1520
1521 }
1522 }
1523 clock_t time1=clock();
1524 cout<<"test_Adj_Complete, time="<<(time1-time0)*1000./CLOCKS_PER_SEC<<"ms"<<endl;
1525}
1526// section begin
1527// section from https://cp-algorithms.com/geometry/tangents-to-two-circles.html#implementation
1528// what it does: calculating tangents to two circles;
1529namespace _MdcGeomSvc_impl{
1530
1531 const double EPS = 1E-9;
1532 /**
1533 * @brief find tangent. should have returned optional<line>, but ok..
1534 *
1535 * @param c
1536 * @param r1
1537 * @param r2
1538 * @param ans
1539 */
1540 void tangents (CLHEP::Hep2Vector c, double r1, double r2, vector<MdcGeomSvc::AdjCandiTgtLine2D> & ans) {
1541 double r = r2 - r1;
1542 double z = squared(c.x()) + squared(c.y());
1543 double d = z - squared(r);
1544 if (d < -EPS) return;
1545 d = sqrt (abs (d));
1546 MdcGeomSvc::AdjCandiTgtLine2D l;
1547 l.a = (c.x() * r + c.y() * d) / z;
1548 l.b = (c.y() * r - c.x() * d) / z;
1549 l.c = r1;
1550 ans.push_back (l);
1551 }
1552
1553 /**
1554 * @brief find all tangents of 2 circles given as point and radius. may return 0,1,2,3,4 lines
1555 *
1556 * @param a
1557 * @param b
1558 * @param ra
1559 * @param rb
1560 * @return vector<MdcGeomSvc::AdjCandiTgtLine2D>
1561 */
1562 vector<MdcGeomSvc::AdjCandiTgtLine2D> tangents (CLHEP::Hep2Vector a, CLHEP::Hep2Vector b, double ra, double rb) {
1563 vector<MdcGeomSvc::AdjCandiTgtLine2D> ans;
1564 for (int i=-1; i<=1; i+=2)
1565 for (int j=-1; j<=1; j+=2)
1566 tangents (b-a, ra*i, rb*j, ans);
1567 for (size_t i=0; i<ans.size(); ++i)
1568 ans[i].c -= ans[i].a * a.x() + ans[i].b * a.y();
1569 return ans;
1570 }
1571}//namespace impl
1572// section end
1573
1574const double MdcGeomSvc::AdjCandiTgtLine2D::INFTY=1./0;
1575/**
1576 * @brief give the lines that is tangent to both drift circles at certain z. require wireA != wireB.
1577 *
1578 * inf is used to rep uninitialized value.
1579 * gives AdjCandiTgtLine2D as the lines. mid of two tangent points are provided with AdjCandiTgtLine2D::location().
1580 * further check AdjCandiTgtLine2D.
1581 *
1582 *
1583 * @param wireA,wireB the 2 wires
1584 * @param driftA,driftB the 2 drifts
1585 * @param expectedZ the z to calc
1586 * @return vector<MdcGeomSvc::AdjCandiTgtLine2D>
1587 */
1588vector<MdcGeomSvc::AdjCandiTgtLine2D> MdcGeomSvc::getAdjacentCandidateWithDriftCircle(const MdcGeoWire* wireA,const MdcGeoWire* wireB, double driftA, double driftB, double expectedZ){
1589
1590 HepPoint3D A1=wireA->Forward();
1591 // HepPoint3D A2(wireA->GetX(expectedZ), wireA->GetY(expectedZ), expectedZ);
1592 HepPoint3D A3=wireA->Backward();
1593 HepPoint3D B1=wireB->Forward(); // is -
1594 // HepPoint3D B2(wireB->GetX(expectedZ), wireB->GetY(expectedZ), expectedZ);
1595 HepPoint3D B3=wireB->Backward();
1596 HepPoint3D directionA=(A3-A1);
1597 HepPoint3D directionB=(B3-B1);
1598 HepPoint3D directionABSurface=directionA.cross(directionB); // perp to the surface.
1599 if (directionABSurface.dot(A1)<0)
1600 directionABSurface=-directionABSurface; // make it outwards
1601 HepPoint3D directionABSurface2D=directionABSurface;
1602 directionABSurface2D.setZ(0);
1603 // double directionABSurface_mag=directionABSurface.mag();
1604 // HepPoint3D directionAlongAB_2D=A2-B2; // along the surface
1605 // double directionAlongAB_2D_mag=directionAlongAB_2D.mag();
1606
1607 // static int runcount=0;
1608 // {
1609 // runcount++;
1610 // if (runcount<5){
1611 // cout<<"debug_out:"<<__LINE__<<" "<<directionAlongAB_2D.x()<<' '<<directionAlongAB_2D.y()<<" "<<directionAlongAB_2D_mag<<endl;
1612 // }
1613 // }
1614 CLHEP::Hep2Vector pA(wireA->GetX(expectedZ),wireA->GetY(expectedZ));
1615 CLHEP::Hep2Vector pB(wireB->GetX(expectedZ),wireB->GetY(expectedZ));
1616 vector<MdcGeomSvc::AdjCandiTgtLine2D> ret=_MdcGeomSvc_impl::tangents(pA,pB, driftA, driftB);
1617 CLHEP::Hep2Vector pMid=(pA+pB)/2;
1618 // line: x0 x + y0 y - x0 ^2 - y0 ^2 =0
1619 // that is, we take the line as the line tangent to the circle passing my point
1620 // double r2=pMid.mag2();
1621 double x0=pMid.x();
1622 double y0=pMid.y();
1623 for (size_t i=0; i<ret.size(); i++) {
1624 double a=ret[i].a;
1625 double b=ret[i].b;
1626 double c=ret[i].c;
1627 // double a2b2=a*a+b*b; //=1
1628 // double mres=(x0*b - y0 * a);
1629 // if (abs(mres)<_MdcGeomSvc_impl::EPS){
1630 // // the root of perp for pMid against line.
1631 // ret[i].x0=AdjCandiTgtLine2D::INFTY;
1632 // ret[i].y0=AdjCandiTgtLine2D::INFTY;
1633 // ret[i].z0=expectedZ;
1634 // continue;
1635 // }
1636 // double x=(r2 * b + c * y0)/mres;
1637 // // double y=-(a*x+c)/b; //what if b ==0.
1638 // double y=(c*x0+a*r2)/(-mres);
1639 double x=(b*b*x0 - a*b*y0 - a*c); // .../a2b2
1640 double y=(-a*b*x0 + a*a*y0 - b*c);
1641
1642 ret[i].x0=x;
1643 ret[i].y0=y;
1644 ret[i].z0=expectedZ;
1645 HepPoint3D directionLine(b,-a,0);// is this right?
1646 if (directionLine.dot(HepPoint3D(x,y,0))<0)
1647 directionLine = -directionLine;
1648 HepPoint3D directionLineRotate90(-directionLine.y(),directionLine.x(),0);// spin pi/2 anti-clock.
1649 //double directionLine_mag=directionLine.mag(); // should be 1;
1650 //ret[i].angle3D=asin(directionLine.dot(directionABSurface )/(directionLine_mag*directionABSurface_mag));
1651 //double cosAngle2D=directionLine.dot(directionAlongAB_2D)/(directionLine_mag*directionAlongAB_2D_mag);
1652
1653 // double sinAngle3D=directionLine.dot(directionABSurface )/(directionABSurface_mag);
1654 // if (sinAngle3D>1 && sinAngle3D < 1 + _MdcGeomSvc_impl::EPS) sinAngle3D=1;
1655 // if (sinAngle3D<-1 && sinAngle3D > -1 - _MdcGeomSvc_impl::EPS)sinAngle3D=-1;
1656 // ret[i].angle3D=asin(sinAngle3D);
1657 ret[i].angle2D=directionLineRotate90.angle(directionABSurface);
1658 ret[i].angle3D=directionLineRotate90.angle(directionABSurface2D);
1659 // double cosAngle2D=directionLine.dot(directionAlongAB_2D)/(directionAlongAB_2D_mag);
1660 // if (cosAngle2D>1 && cosAngle2D < 1+ _MdcGeomSvc_impl::EPS)cosAngle2D=1;
1661 // if (cosAngle2D<-1 && cosAngle2D > -1 - _MdcGeomSvc_impl::EPS) cosAngle2D=-1;
1662 // // in test, rounding errors can make the cosAngle>1 when it is supposed to be 1. we revert that effect here.
1663 // ret[i].angle2D=acos(cosAngle2D);
1664 // {
1665 // if (runcount<5){
1666 // cout<<"debug_out:"<<__LINE__<<" "<<directionLine.x()<<' '<<directionLine.y()<<" "<<directionLine_mag<< " "<<directionLine.dot(directionAlongAB_2D)/(directionLine_mag*directionAlongAB_2D_mag)<<endl;
1667 // }
1668 // }
1669
1670 }// adding point of Hit.
1671
1672 return ret;
1673
1674}
1675
1676
1677/**
1678 * @brief give the lines that is tangent to both drift circles at the approximate z that the wires cross. require wireA != wireB.
1679 *
1680 * calls getAdjacentCandidateWithDriftCircle() with given Z. z calculated from center of bounds given in getAdjacentIntersectionBounds();
1681 * if no bound avaliable (the 2 wires has no neighbor), no line is returnd
1682 *
1683 * @param wireA,wireB the 2 wires
1684 * @param driftA,driftB the 2 drifts
1685 * @return vector<CLHEP::Hep2Vector>
1686 */
1687vector<MdcGeomSvc::AdjCandiTgtLine2D> MdcGeomSvc::getAdjacentCandidateWithDriftCircle(const MdcGeoWire* wireA,const MdcGeoWire* wireB, double driftA, double driftB){
1688 std::vector<double> bound=getAdjacentIntersectionBounds(wireA, wireB);
1689 //[phiL,phiU, zL,zU]
1690 double expectedZ;
1691 if (bound.size()==4)expectedZ=(bound[3]+bound[2])/2;
1692 else {
1693 vector<MdcGeomSvc::AdjCandiTgtLine2D> ret;
1694 return ret ;
1695 }
1696 return getAdjacentCandidateWithDriftCircle(wireA, wireB, driftA, driftB,expectedZ);
1697}
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
nPhi
const int nPhi
Definition BesTofDigitizerEcV3.hh:37
Phi
double Phi(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:99
length
double length
Definition CgemLineFit.cxx:121
n
const Int_t n
Definition DataBase/tau_mode.c:65
phi2
Double_t phi2
Definition DataBase/tau_mode.c:7
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
phi1
Double_t phi1
Definition DataBase/tau_mode.c:7
wireNo
const int wireNo
Definition DedxCalibWireGain.cxx:18
abs
double abs(const EvtComplex &c)
Definition EvtComplex.hh:212
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition Gam4pikp.cxx:37
DBL_MAX
#define DBL_MAX
Definition KalFitAlg.h:13
M_PI
#define M_PI
Definition TConstant.h:4
MdcGeoWire::GetX
double GetX(const double z) const
Definition MdcGeomSvc-00-01-42/src/MdcGeoWire.cxx:49
MdcGeoWire::GetY
double GetY(const double z) const
Definition MdcGeomSvc-00-01-42/src/MdcGeoWire.cxx:53
MdcGeoWire::GetNeighborType4
vector< MdcGeoWire * > GetNeighborType4(void) const
Definition MdcGeomSvc-00-01-42/src/MdcGeoWire.cxx:87
MdcGeoWire::WirePhi
double WirePhi(double z) const
Definition MdcGeomSvc-00-01-42/src/MdcGeoWire.cxx:146
MdcGeomSvc::MdcGeomSvc
MdcGeomSvc(const std::string &name, ISvcLocator *svcloc)
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:30
MdcGeomSvc::getAdjacentCandidateWithDriftCircle
vector< AdjCandiTgtLine2D > getAdjacentCandidateWithDriftCircle(const MdcGeoWire *wireA, const MdcGeoWire *wireB, double driftA, double driftB)
give the lines that is tangent to both drift circles at the approximate z that the wires cross....
Definition MdcGeomSvc-00-01-42/src/MdcGeomSvc.cxx:1687
MdcGeomSvc::SuperLayer
const MdcGeoSuper *const SuperLayer(unsigned id)
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:794
MdcGeomSvc::Wire
const MdcGeoWire *const Wire(unsigned id)
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:770
MdcGeomSvc::GeneralLayer
const MdcGeoGeneral *const GeneralLayer(unsigned id)
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:802
MdcGeomSvc::getAdjacentIntersectionPointFast
HepPoint3D getAdjacentIntersectionPointFast(const MdcGeoWire *wireA, const MdcGeoWire *wireB)
find the z where points on wireA, wireB has same phi. point is given as avg(A(z),B(z))
Definition MdcGeomSvc-00-01-42/src/MdcGeomSvc.cxx:1249
MdcGeomSvc::handle
void handle(const Incident &inc)
this handle function is prepared for special use
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:743
MdcGeomSvc::getAdjacentIntersectionPointSlower
HepPoint3D getAdjacentIntersectionPointSlower(const MdcGeoWire *wireA, const MdcGeoWire *wireB)
Definition MdcGeomSvc-00-01-42/src/MdcGeomSvc.cxx:1275
MdcGeomSvc::getAdjacentIntersectionBounds
vector< double > getAdjacentIntersectionBounds(const MdcGeoWire *wireA, const MdcGeoWire *wireB)
get bounding box in \varphi and z of touching shape of 2 adjacent layer wires.
Definition MdcGeomSvc-00-01-42/src/MdcGeomSvc.cxx:1298
MdcGeomSvc::queryInterface
virtual StatusCode queryInterface(const InterfaceID &riid, void **ppvUnknown)
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:63
MdcGeomSvc::End
const MdcGeoEnd *const End(unsigned id)
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:815
MdcGeomSvc::Layer
const MdcGeoLayer *const Layer(unsigned id)
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:786
MdcGeomSvc::getAdjacentIntersectionShape
vector< HepPoint3D > getAdjacentIntersectionShape(const MdcGeoWire *wireA, const MdcGeoWire *wireB)
get touching shape of 2 adjacent layer wires. return vertexs of intersection shape if the wires' se...
Definition MdcGeomSvc-00-01-42/src/MdcGeomSvc.cxx:1045
MdcGeomSvc::Misc
const MdcGeoMisc *const Misc(void)
Definition bak_MdcGeomSvc-00-01-39/src/MdcGeomSvc.cxx:810
_MdcGeomSvc_impl
Definition MdcGeomSvc-00-01-42/src/MdcGeomSvc.cxx:1529
_MdcGeomSvc_impl::EPS
const double EPS
Definition MdcGeomSvc-00-01-42/src/MdcGeomSvc.cxx:1531
_MdcGeomSvc_impl::tangents
void tangents(CLHEP::Hep2Vector c, double r1, double r2, vector< MdcGeomSvc::AdjCandiTgtLine2D > &ans)
find tangent. should have returned optional<line>, but ok..
Definition MdcGeomSvc-00-01-42/src/MdcGeomSvc.cxx:1540
rb
Definition TConstant.h:3
rb::R
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition TUtil.h:27
MdcGeomSvc::AdjCandiTgtLine2D
forms a line in x-y space, ax + by + c = 0
Definition MdcGeomSvc-00-01-42/MdcGeomSvc/MdcGeomSvc.h:58
pi
const float pi
Definition vector3.h:133