1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //*************************************************************************
18 // Class for flat cables
20 // Ludovic Gaudichet gaudichet@to.infn.it
21 //*************************************************************************
25 // General Root includes
26 //#include <Riostream.h>
30 // Root Geometry includes
31 #include <TGeoManager.h>
32 #include <TGeoVolume.h>
35 #include <TGeoMatrix.h>
38 #include "AliITSv11GeomCableFlat.h"
41 ClassImp(AliITSv11GeomCableFlat)
43 //________________________________________________________________________
44 AliITSv11GeomCableFlat::AliITSv11GeomCableFlat():
51 for (Int_t i=0; i<fgkCableMaxLayer ; i++) {
59 //________________________________________________________________________
60 AliITSv11GeomCableFlat::
61 AliITSv11GeomCableFlat(const char* name, Double_t width, Double_t thick) :
62 AliITSv11GeomCable(name),
67 // standard constructor
68 for (Int_t i=0; i<fgkCableMaxLayer ; i++) {
76 //________________________________________________________________________
77 AliITSv11GeomCableFlat::AliITSv11GeomCableFlat(const AliITSv11GeomCableFlat &s) :
78 AliITSv11GeomCable(s),fWidth(s.fWidth),fThick(s.fThick),fNlayer(s.fNlayer)
81 for (Int_t i=0; i<s.fNlayer; i++) {
82 fLayThickness[i] = s.fLayThickness[i];
83 fTranslation[i] = s.fTranslation[i];
84 fLayMedia[i] = s.fLayMedia[i];
85 fLayColor[i] = s.fLayColor[i];
89 //________________________________________________________________________
90 AliITSv11GeomCableFlat& AliITSv11GeomCableFlat::
91 operator=(const AliITSv11GeomCableFlat &s) {
92 // Assignment operator
93 // Not fully inplemented yet !!!
95 if(&s == this) return *this;
100 for (Int_t i=0; i<s.fNlayer; i++) {
101 fLayThickness[i] = s.fLayThickness[i];
102 fTranslation[i] = s.fTranslation[i];
103 fLayMedia[i] = s.fLayMedia[i];
104 fLayColor[i] = s.fLayColor[i];
109 //________________________________________________________________________
110 Int_t AliITSv11GeomCableFlat::GetPoint( Int_t iCheckPt, Double_t *coord)
112 // Get the correct point #iCheckPt
113 TVectorD *coordVector =(TVectorD *)fPointArray.At(2*iCheckPt);
115 #if ROOT_VERSION_CODE < ROOT_VERSION(4,0,0)
116 CopyFrom(coord, coordVector->GetElements());
118 CopyFrom(coord, coordVector->GetMatrixArray());
127 //________________________________________________________________________
128 Int_t AliITSv11GeomCableFlat::GetVect( Int_t iCheckPt, Double_t *coord)
130 // Get the correct vect corresponding to point #iCheckPt
132 TVectorD *coordVector =(TVectorD *)fPointArray.At(2*iCheckPt+1);
134 #if ROOT_VERSION_CODE < ROOT_VERSION(4,0,0)
135 CopyFrom(coord, coordVector->GetElements());
137 CopyFrom(coord, coordVector->GetMatrixArray());
146 //________________________________________________________________________
147 void AliITSv11GeomCableFlat::AddCheckPoint( TGeoVolume *vol, Int_t iCheckPt,
148 Double_t *coord, Double_t *orthVect)
151 // Add a check point. In the fPointArray, the point is at i and its vector
155 // if (iCheckPt>=fVolumeArray.GetEntriesFast()) {
156 // fVolumeArray.AddLast(vol);
157 // TVectorD *point = new TVectorD(3,coord);
158 // TVectorD *vect = new TVectorD(3,orthVect);
159 // fPointArray.AddLast(point);
160 // fPointArray.AddLast(vect);
162 // } else if ((iCheckPt >= 0)&&(iCheckPt < fVolumeArray.GetEntriesFast())) {
163 // fVolumeArray.AddAt(vol, iCheckPt);
164 // TVectorD *point = new TVectorD(3,coord);
165 // TVectorD *vect = new TVectorD(3,orthVect);
166 // fPointArray.AddAt(point, iCheckPt*2 );
167 // fPointArray.AddAt(vect, iCheckPt*2+1);
169 fVolumeArray.AddAtAndExpand(vol, iCheckPt);
170 TVectorD *point = new TVectorD(3,coord);
171 TVectorD *vect = new TVectorD(3,orthVect);
172 fPointArray.AddAtAndExpand(point, iCheckPt*2 );
173 fPointArray.AddAtAndExpand(vect, iCheckPt*2+1);
176 //________________________________________________________________________
177 void AliITSv11GeomCableFlat::PrintCheckPoints() const {
178 // print all check points of the cable
179 printf(" ---\n Printing all check points of the flat cable\n");
180 for (Int_t i = 0; i<fVolumeArray.GetEntriesFast(); i++) {
182 if (GetPoint( i, coord))
183 printf(" ( %.2f, %.2f, %.2f )\n", coord[0], coord[1], coord[2]);
187 //________________________________________________________________________
188 TGeoVolume* AliITSv11GeomCableFlat::CreateAndInsertCableSegment(Int_t p2,
192 // Creates a cable segment between points p1 and p2.
193 // Rotation is the eventual rotation of the flat cable
194 // along its length axis
196 // The segment volume is created inside the volume containing point2
197 // Therefore this segment should be defined in this volume only.
198 // I mean here that, if the previous point is in another volume,
199 // it should be just at the border between the 2 volumes. Also the
200 // orientation vector of the previous point should be orthogonal to
201 // the surface between the 2 volumes.
204 if (fInitialNode==0) {
205 TObjArray *nodes = gGeoManager->GetListOfNodes();
206 if (nodes->GetEntriesFast()==0) return 0;
207 mainNode = (TGeoNode *) nodes->UncheckedAt(0);
209 mainNode = fInitialNode;
213 TGeoVolume *p2Vol = GetVolume(p2);
214 TGeoVolume *p1Vol = GetVolume(p1);
216 ResetCheckDaughter();
218 if (! CheckDaughter(mainNode)) {
219 printf("Error::volume containing point is not visible in node tree!\n");
223 Double_t coord1[3], coord2[3], vect1[3], vect2[3];
224 //=================================================
225 // Get p1 position in the systeme of p2
228 Int_t p1nodeInd[fgkCableMaxNodeLevel];
229 for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p1nodeInd[i]=fNodeInd[i];
230 Int_t p1volLevel = 0;
231 while (p1nodeInd[p1volLevel]!=-1) p1volLevel++;
234 ResetCheckDaughter();
236 if (! CheckDaughter(mainNode)) {
237 printf("Error::volume containing point is not visible in node tree!\n");
240 Int_t p2nodeInd[fgkCableMaxNodeLevel];
241 for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p2nodeInd[i]=fNodeInd[i];
242 Int_t commonMotherLevel = 0;
243 while (p1nodeInd[commonMotherLevel]==fNodeInd[commonMotherLevel])
246 Int_t p2volLevel = 0;
247 while (fNodeInd[p2volLevel]!=-1) p2volLevel++;
250 // Get coord and vect of p1 in the common mother reference system
251 if (! GetCheckPoint(p1, 0, p1volLevel-commonMotherLevel, coord1) )
253 if (! GetCheckVect( p1, 0, p1volLevel-commonMotherLevel, vect1) )
256 // Translate them in the reference system of the volume containing p2
257 TGeoNode *pathNode[fgkCableMaxNodeLevel];
258 pathNode[0] = mainNode;
259 for (Int_t i=0; i<=p2volLevel; i++) {
260 pathNode[i+1] = pathNode[i]->GetDaughter(p2nodeInd[i]);
262 Double_t globalCoord1[3] = {coord1[0], coord1[1], coord1[2]};
263 Double_t globalVect1[3] = {vect1[0], vect1[1], vect1[2]};
265 for (Int_t i = commonMotherLevel+1; i <= p2volLevel; i++) {
266 pathNode[i+1]->GetMatrix()->MasterToLocal(globalCoord1, coord1);
267 pathNode[i+1]->GetMatrix()->MasterToLocalVect(globalVect1, vect1);
268 CopyFrom(globalCoord1, coord1);
269 CopyFrom(globalVect1, vect1);
272 if (! GetCheckPoint(p1, 0, 0, coord1) ) return 0;
273 if (! GetCheckVect(p1, 0, 0, vect1) ) return 0;
276 //=================================================
277 // Get p2 position in the systeme of p2
278 if (! GetCheckPoint(p2, 0, 0, coord2) ) return 0;
279 if (! GetCheckVect(p2, 0, 0, vect2) ) return 0;
281 Double_t cx = (coord1[0]+coord2[0])/2;
282 Double_t cy = (coord1[1]+coord2[1])/2;
283 Double_t cz = (coord1[2]+coord2[2])/2;
284 Double_t dx = coord2[0]-coord1[0];
285 Double_t dy = coord2[1]-coord1[1];
286 Double_t dz = coord2[2]-coord1[2];
288 //=================================================
289 // Positionning of the segment between the 2 points
290 if (TMath::Abs(dy)<1e-231) dy = 1e-231;
291 if (TMath::Abs(dz)<1e-231) dz = 1e-231;
292 //Double_t angleRot1 = -TMath::ATan(dx/dy);
293 //Double_t planDiagL = -TMath::Sqrt(dy*dy+dx*dx);
294 //if (dy<0) planDiagL = -planDiagL;
295 //Double_t angleRotDiag = TMath::ATan(planDiagL/dz);
297 Double_t angleRot1 = -TMath::ATan2(dx,dy);
298 Double_t planDiagL = TMath::Sqrt(dy*dy+dx*dx);
299 Double_t angleRotDiag = -TMath::ATan2(planDiagL,dz);
300 //--- (Calculate rotation of segment on the Z axis)
301 //-- Here I'm trying to calculate the rotation to be applied in
302 //-- order to match as closer as possible this segment and the
304 //-- It seems that some times it doesn't work ...
305 Double_t angleRotZ = 0;
306 TGeoRotation rotTemp("",angleRot1*TMath::RadToDeg(),
307 angleRotDiag*TMath::RadToDeg(), rotation);
308 Double_t localX[3] = {0,1,0};
310 rotTemp.LocalToMasterVect(localX, globalX);
311 CopyFrom(localX, globalX);
312 GetCheckVect(localX, p2Vol, 0, fgkCableMaxNodeLevel+1, globalX);
313 Double_t orthVect[3];
314 GetCheckVect(vect1, p2Vol, 0, fgkCableMaxNodeLevel+1, orthVect);
316 Double_t orthVectNorm2 = ScalProd(orthVect,orthVect);
317 Double_t alpha1 = ScalProd(fPreviousX,orthVect)/orthVectNorm2;
318 Double_t alpha2 = ScalProd(globalX,orthVect)/orthVectNorm2;
319 Double_t globalX1p[3], globalX2p[3];
320 globalX1p[0] = fPreviousX[0] - alpha1*orthVect[0];
321 globalX1p[1] = fPreviousX[1] - alpha1*orthVect[1];
322 globalX1p[2] = fPreviousX[2] - alpha1*orthVect[2];
323 globalX2p[0] = globalX[0] - alpha2*orthVect[0];
324 globalX2p[1] = globalX[1] - alpha2*orthVect[1];
325 globalX2p[2] = globalX[2] - alpha2*orthVect[2];
326 //-- now I'm searching the 3th vect which makes an orthogonal base
327 //-- with orthVect and globalX1p ...
328 Double_t nulVect[3] = {0,0,0};
330 TMath::Normal2Plane(nulVect, orthVect, globalX1p, axis3);
331 Double_t globalX1pNorm2 = ScalProd(globalX1p, globalX1p);
332 Double_t beta = ScalProd(globalX2p, globalX1p)/globalX1pNorm2;
333 Double_t gamma = ScalProd(globalX2p, axis3);
334 angleRotZ = (TMath::ATan2(1,0) - TMath::ATan2(beta, gamma))
337 // cout << "!!!!!!!!!!!!!!!!!!! angle = " <<angleRotZ << endl;
338 CopyFrom(fPreviousX, globalX);
340 Double_t localVect1[3], localVect2[3];
341 TGeoRotation rot("",angleRot1*TMath::RadToDeg(),
342 angleRotDiag*TMath::RadToDeg(),
344 // rotation-angleRotZ);
345 // since angleRotZ doesn't always work, I won't use it ...
347 rot.MasterToLocalVect(vect1, localVect1);
348 rot.MasterToLocalVect(vect2, localVect2);
350 //=================================================
351 // Create the segment and add it to the mother volume
352 TGeoVolume *vCableSegB = CreateSegment(coord1, coord2,
353 localVect1, localVect2);
354 // TGeoVolume *vCableSegB = CreateBoxSegment(coord1, coord2);
356 TGeoRotation rotArbSeg("", 0, 90, 0);
357 rotArbSeg.MultiplyBy(&rot, kFALSE);
358 TGeoTranslation trans("",cx, cy, cz);
359 TGeoCombiTrans *combiB = new TGeoCombiTrans(trans, rotArbSeg);
360 p2Vol->AddNode(vCableSegB, p2, combiB);
361 //=================================================;
364 printf("---\n Cable segment points : ");
365 printf("%f, %f, %f\n",coord1[0], coord1[1], coord1[2]);
366 printf("%f, %f, %f\n",coord2[0], coord2[1], coord2[2]);
369 // #include <TGeoSphere.h>
370 // TGeoMedium *airSDD = gGeoManager->GetMedium("ITS_AIR$");
371 // TGeoSphere *sphere = new TGeoSphere(0, 0.05);
372 // TGeoVolume *vSphere = new TGeoVolume("", sphere, airSDD);
373 // TGeoTranslation *trC = new TGeoTranslation("", cx, cy, cz);
374 // TGeoTranslation *tr1 = new TGeoTranslation("",coord1[0],
375 // coord1[1],coord1[2]);
376 // TGeoTranslation *tr2 = new TGeoTranslation("",coord2[0],
377 // coord2[1],coord2[2]);
378 // p2Vol->AddNode(vSphere, p2*3-2, trC);
379 // p2Vol->AddNode(vSphere, p2*3-1, tr1);
380 // p2Vol->AddNode(vSphere, p2*3 , tr2);
381 if (ct) *ct = combiB;
385 //________________________________________________________________________
386 TGeoVolume* AliITSv11GeomCableFlat::CreateAndInsertBoxCableSegment(Int_t p2,
390 // This function is to be use only when the segment has the shape
391 // of a simple box, i.e. the normal vector to its end is perpendicular
392 // to the segment own axis
393 // Creates a cable segment between points p1 and p2.
394 // Rotation is the eventual rotation of the flat cable
395 // along its length axis
397 // The segment volume is created inside the volume containing point2
398 // Therefore this segment should be defined in this volume only.
399 // I mean here that, if the previous point is in another volume,
400 // it should be just at the border between the 2 volumes. Also the
401 // orientation vector of the previous point should be orthogonal to
402 // the surface between the 2 volumes.
405 if (fInitialNode==0) {
406 TObjArray *nodes = gGeoManager->GetListOfNodes();
407 if (nodes->GetEntriesFast()==0) return 0;
408 mainNode = (TGeoNode *) nodes->UncheckedAt(0);
410 mainNode = fInitialNode;
414 TGeoVolume *p2Vol = GetVolume(p2);
415 TGeoVolume *p1Vol = GetVolume(p1);
417 ResetCheckDaughter();
419 if (! CheckDaughter(mainNode)) {
420 printf("Error::volume containing point is not visible in node tree!\n");
424 Double_t coord1[3], coord2[3], vect1[3], vect2[3];
425 //=================================================
426 // Get p1 position in the systeme of p2
429 Int_t p1nodeInd[fgkCableMaxNodeLevel];
430 for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p1nodeInd[i]=fNodeInd[i];
431 Int_t p1volLevel = 0;
432 while (p1nodeInd[p1volLevel]!=-1) p1volLevel++;
435 ResetCheckDaughter();
437 if (! CheckDaughter(mainNode)) {
438 printf("Error::volume containing point is not visible in node tree!\n");
441 Int_t p2nodeInd[fgkCableMaxNodeLevel];
442 for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p2nodeInd[i]=fNodeInd[i];
443 Int_t commonMotherLevel = 0;
444 while (p1nodeInd[commonMotherLevel]==fNodeInd[commonMotherLevel])
447 Int_t p2volLevel = 0;
448 while (fNodeInd[p2volLevel]!=-1) p2volLevel++;
451 // Get coord and vect of p1 in the common mother reference system
452 if (! GetCheckPoint(p1, 0, p1volLevel-commonMotherLevel, coord1) )
454 if (! GetCheckVect( p1, 0, p1volLevel-commonMotherLevel, vect1) )
457 // Translate them in the reference system of the volume containing p2
458 TGeoNode *pathNode[fgkCableMaxNodeLevel];
459 pathNode[0] = mainNode;
460 for (Int_t i=0; i<=p2volLevel; i++) {
461 pathNode[i+1] = pathNode[i]->GetDaughter(p2nodeInd[i]);
463 Double_t globalCoord1[3] = {coord1[0], coord1[1], coord1[2]};
464 Double_t globalVect1[3] = {vect1[0], vect1[1], vect1[2]};
466 for (Int_t i = commonMotherLevel+1; i <= p2volLevel; i++) {
467 pathNode[i+1]->GetMatrix()->MasterToLocal(globalCoord1, coord1);
468 pathNode[i+1]->GetMatrix()->MasterToLocalVect(globalVect1, vect1);
469 CopyFrom(globalCoord1, coord1);
470 CopyFrom(globalVect1, vect1);
473 if (! GetCheckPoint(p1, 0, 0, coord1) ) return 0;
474 if (! GetCheckVect(p1, 0, 0, vect1) ) return 0;
477 //=================================================
478 // Get p2 position in the systeme of p2
479 if (! GetCheckPoint(p2, 0, 0, coord2) ) return 0;
480 if (! GetCheckVect(p2, 0, 0, vect2) ) return 0;
482 Double_t cx = (coord1[0]+coord2[0])/2;
483 Double_t cy = (coord1[1]+coord2[1])/2;
484 Double_t cz = (coord1[2]+coord2[2])/2;
485 Double_t dx = coord2[0]-coord1[0];
486 Double_t dy = coord2[1]-coord1[1];
487 Double_t dz = coord2[2]-coord1[2];
489 //=================================================
490 // Positionning of the segment between the 2 points
491 if (TMath::Abs(dy)<1e-231) dy = 1e-231;
492 if (TMath::Abs(dz)<1e-231) dz = 1e-231;
493 //Double_t angleRot1 = -TMath::ATan(dx/dy);
494 //Double_t planDiagL = -TMath::Sqrt(dy*dy+dx*dx);
495 //if (dy<0) planDiagL = -planDiagL;
496 //Double_t angleRotDiag = TMath::ATan(planDiagL/dz);
498 Double_t angleRot1 = -TMath::ATan2(dx,dy);
499 Double_t planDiagL = TMath::Sqrt(dy*dy+dx*dx);
500 Double_t angleRotDiag = -TMath::ATan2(planDiagL,dz);
501 //--- (Calculate rotation of segment on the Z axis)
502 //-- Here I'm trying to calculate the rotation to be applied in
503 //-- order to match as closer as possible this segment and the
505 //-- It seems that some times it doesn't work ...
506 Double_t angleRotZ = 0;
507 TGeoRotation rotTemp("",angleRot1*TMath::RadToDeg(),
508 angleRotDiag*TMath::RadToDeg(), rotation);
509 Double_t localX[3] = {0,1,0};
511 rotTemp.LocalToMasterVect(localX, globalX);
512 CopyFrom(localX, globalX);
513 GetCheckVect(localX, p2Vol, 0, fgkCableMaxNodeLevel+1, globalX);
514 Double_t orthVect[3];
515 GetCheckVect(vect1, p2Vol, 0, fgkCableMaxNodeLevel+1, orthVect);
517 Double_t orthVectNorm2 = ScalProd(orthVect,orthVect);
518 Double_t alpha1 = ScalProd(fPreviousX,orthVect)/orthVectNorm2;
519 Double_t alpha2 = ScalProd(globalX,orthVect)/orthVectNorm2;
520 Double_t globalX1p[3], globalX2p[3];
521 globalX1p[0] = fPreviousX[0] - alpha1*orthVect[0];
522 globalX1p[1] = fPreviousX[1] - alpha1*orthVect[1];
523 globalX1p[2] = fPreviousX[2] - alpha1*orthVect[2];
524 globalX2p[0] = globalX[0] - alpha2*orthVect[0];
525 globalX2p[1] = globalX[1] - alpha2*orthVect[1];
526 globalX2p[2] = globalX[2] - alpha2*orthVect[2];
527 //-- now I'm searching the 3th vect which makes an orthogonal base
528 //-- with orthVect and globalX1p ...
529 Double_t nulVect[3] = {0,0,0};
531 TMath::Normal2Plane(nulVect, orthVect, globalX1p, axis3);
532 Double_t globalX1pNorm2 = ScalProd(globalX1p, globalX1p);
533 Double_t beta = ScalProd(globalX2p, globalX1p)/globalX1pNorm2;
534 Double_t gamma = ScalProd(globalX2p, axis3);
535 angleRotZ = (TMath::ATan2(1,0) - TMath::ATan2(beta, gamma))
538 // cout << "!!!!!!!!!!!!!!!!!!! angle = " <<angleRotZ << endl;
539 CopyFrom(fPreviousX, globalX);
541 Double_t localVect1[3], localVect2[3];
542 TGeoRotation rot("",angleRot1*TMath::RadToDeg(),
543 angleRotDiag*TMath::RadToDeg(),
545 // rotation-angleRotZ);
546 // since angleRotZ doesn't always work, I won't use it ...
548 rot.MasterToLocalVect(vect1, localVect1);
549 rot.MasterToLocalVect(vect2, localVect2);
551 //=================================================
552 // Create the segment and add it to the mother volume
553 TGeoVolume *vCableSegB = CreateBoxSegment(coord1, coord2);
555 TGeoRotation rotArbSeg("", 0, 90, 0);
556 rotArbSeg.MultiplyBy(&rot, kFALSE);
557 TGeoTranslation trans("",cx, cy, cz);
558 TGeoCombiTrans *combiB = new TGeoCombiTrans(trans, rotArbSeg);
559 p2Vol->AddNode(vCableSegB, p2, combiB);
560 //=================================================;
563 printf("---\n Cable segment points : ");
564 printf("%f, %f, %f\n",coord1[0], coord1[1], coord1[2]);
565 printf("%f, %f, %f\n",coord2[0], coord2[1], coord2[2]);
568 if (ct) *ct = combiB;
572 //________________________________________________________________________
573 TGeoVolume* AliITSv11GeomCableFlat::CreateAndInsertCableCylSegment(Int_t p2,
577 // Create a flat cable segment with a curvature between points p1 and p2.
578 // The radius and position of the curve is defined by the
579 // perpendicular vector of point p2 (the orientation of this vector
580 // and the position of the 2 check points are enough to completely
582 // Rotation is the eventual rotation of the flat cable
583 // along its length axis
587 if (fInitialNode==0) {
588 TObjArray *nodes = gGeoManager->GetListOfNodes();
589 if (nodes->GetEntriesFast()==0) return 0;
590 mainNode = (TGeoNode *) nodes->UncheckedAt(0);
592 mainNode = fInitialNode;
596 TGeoVolume *p1Vol = GetVolume(p1);
597 TGeoVolume *p2Vol = GetVolume(p2);
599 ResetCheckDaughter();
601 if (! CheckDaughter(mainNode)) {
602 printf("Error::volume containing point is not visible in node tree!\n");
606 Double_t coord1[3], coord2[3], vect1[3], vect2[3];
607 //=================================================
608 // Get p1 position in the systeme of p2
611 Int_t p1nodeInd[fgkCableMaxNodeLevel];
612 for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p1nodeInd[i]=fNodeInd[i];
613 Int_t p1volLevel = 0;
614 while (p1nodeInd[p1volLevel]!=-1) p1volLevel++;
617 ResetCheckDaughter();
619 if (! CheckDaughter(mainNode)) {
620 printf("Error::volume containing point is not visible in node tree!\n");
623 Int_t p2nodeInd[fgkCableMaxNodeLevel];
624 for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p2nodeInd[i]=fNodeInd[i];
625 Int_t commonMotherLevel = 0;
626 while (p1nodeInd[commonMotherLevel]==fNodeInd[commonMotherLevel])
629 Int_t p2volLevel = 0;
630 while (fNodeInd[p2volLevel]!=-1) p2volLevel++;
633 // Get coord and vect of p1 in the common mother reference system
634 GetCheckPoint(p1, 0, p1volLevel-commonMotherLevel, coord1);
635 GetCheckVect( p1, 0, p1volLevel-commonMotherLevel, vect1);
636 // Translate them in the reference system of the volume containing p2
637 TGeoNode *pathNode[fgkCableMaxNodeLevel];
638 pathNode[0] = mainNode;
639 for (Int_t i=0; i<=p2volLevel; i++) {
640 pathNode[i+1] = pathNode[i]->GetDaughter(p2nodeInd[i]);
642 Double_t globalCoord1[3] = {coord1[0], coord1[1], coord1[2]};
643 Double_t globalVect1[3] = {vect1[0], vect1[1], vect1[2]};
645 for (Int_t i = commonMotherLevel+1; i<=p2volLevel; i++) {
646 pathNode[i+1]->GetMatrix()->MasterToLocal(globalCoord1, coord1);
647 pathNode[i+1]->GetMatrix()->MasterToLocalVect(globalVect1, vect1);
648 CopyFrom(globalCoord1, coord1);
649 CopyFrom(globalVect1, vect1);
652 GetCheckPoint(p1, 0, 0, coord1);
653 GetCheckVect(p1, 0, 0, vect1);
656 //=================================================
657 // Get p2 position in the systeme of p2
658 GetCheckPoint(p2, 0, 0, coord2);
659 GetCheckVect(p2, 0, 0, vect2);
661 Double_t cx = (coord1[0]+coord2[0])/2;
662 Double_t cy = (coord1[1]+coord2[1])/2;
663 Double_t cz = (coord1[2]+coord2[2])/2;
664 Double_t dx = coord2[0]-coord1[0];
665 Double_t dy = coord2[1]-coord1[1];
666 Double_t dz = coord2[2]-coord1[2];
667 Double_t length = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
669 //=================================================
670 // Positionning of the segment between the 2 points
671 if ((dy<1e-31)&&(dy>0)) dy = 1e-31;
672 if ((dz<1e-31)&&(dz>0)) dz = 1e-31;
673 if ((dy>-1e-31)&&(dy<0)) dy = -1e-31;
674 if ((dz>-1e-31)&&(dz<0)) dz = -1e-31;
676 Double_t angleRot1 = -TMath::ATan2(dx,dy);
677 Double_t planDiagL = TMath::Sqrt(dy*dy+dx*dx);
678 Double_t angleRotDiag = -TMath::ATan2(planDiagL,dz);
680 TGeoRotation rotTorusTemp("",angleRot1*TMath::RadToDeg(),
681 angleRotDiag*TMath::RadToDeg(),0);
682 TGeoRotation rotTorusToZ("",0,90,0);
683 rotTorusTemp.MultiplyBy(&rotTorusToZ, kTRUE);
684 Double_t localVect2[3];
685 rotTorusTemp.MasterToLocalVect(vect2, localVect2);
686 if (localVect2[1]<0) {
687 localVect2[0] = -localVect2[0];
688 localVect2[1] = -localVect2[1];
689 localVect2[2] = -localVect2[2];
691 Double_t normVect2 = TMath::Sqrt(localVect2[0]*localVect2[0]+
692 localVect2[1]*localVect2[1]+
693 localVect2[2]*localVect2[2]);
694 Double_t axisX[3] = {1,0,0};
695 Double_t cosangleTorusSeg = (localVect2[0]*axisX[0]+
696 localVect2[1]*axisX[1]+
697 localVect2[2]*axisX[2])/normVect2;
698 Double_t angleTorusSeg = TMath::ACos(cosangleTorusSeg)*TMath::RadToDeg();
699 TGeoRotation rotTorus("",angleRot1*TMath::RadToDeg(),
700 angleRotDiag*TMath::RadToDeg(),
701 45-angleTorusSeg+rotation);
702 //180-angleTorusSeg+rotation);
703 rotTorus.MultiplyBy(&rotTorusToZ, kTRUE);
704 rotTorus.MasterToLocalVect(vect2, localVect2);
705 if (localVect2[1]<0) {
706 localVect2[0] = -localVect2[0];
707 localVect2[1] = -localVect2[1];
708 localVect2[2] = -localVect2[2];
710 normVect2 = TMath::Sqrt(localVect2[0]*localVect2[0]+
711 localVect2[1]*localVect2[1]+
712 localVect2[2]*localVect2[2]);
713 Double_t axisY[3] = {0,1,0};
714 Double_t cosPhi = (localVect2[0]*axisY[0]+localVect2[1]*axisY[1]+
715 localVect2[2]*axisY[2])/normVect2;
716 Double_t torusPhi1 = TMath::ACos(cosPhi);
717 Double_t torusR = (length/2)/TMath::Sin(torusPhi1);
718 torusPhi1 = torusPhi1*TMath::RadToDeg();
719 Double_t perpLength = TMath::Sqrt(torusR*torusR-length*length/4);
720 Double_t localTransT[3] = {-perpLength,0,0};
721 Double_t globalTransT[3];
722 rotTorus.LocalToMasterVect(localTransT, globalTransT);
723 TGeoTranslation transTorus("",cx+globalTransT[0],cy+globalTransT[1],
726 TGeoCombiTrans *combiTorus = new TGeoCombiTrans(transTorus, rotTorus);
728 //=================================================
729 // Create the segment and add it to the mother volume
730 TGeoVolume *vCableSegT = CreateCylSegment(torusPhi1, torusR);
731 p2Vol->AddNode(vCableSegT, p2, combiTorus);
734 printf("---\n Cable segment points : ");
735 printf("%f, %f, %f\n",coord1[0], coord1[1], coord1[2]);
736 printf("%f, %f, %f\n",coord2[0], coord2[1], coord2[2]);
739 if (ct) *ct = combiTorus;
744 //________________________________________________________________________
745 TGeoVolume *AliITSv11GeomCableFlat::CreateSegment( Double_t *coord1,
747 Double_t *localVect1,
748 Double_t *localVect2 )
751 //=================================================
752 // Calculate segment "deformation"
753 Double_t dx = coord2[0]-coord1[0];
754 Double_t dy = coord2[1]-coord1[1];
755 Double_t dz = coord2[2]-coord1[2];
756 Double_t length = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
758 Double_t cosTheta1 = -1./TMath::Sqrt( 1 + localVect1[0]*localVect1[0]
759 /localVect1[2]/localVect1[2] );
760 Double_t cosTheta2 = 1./TMath::Sqrt( 1 + localVect2[0]*localVect2[0]
761 /localVect2[2]/localVect2[2] );
762 if (localVect1[2]<0) cosTheta1 = -cosTheta1;
763 if (localVect2[2]<0) cosTheta2 = -cosTheta2;
765 Double_t dL1 = 0.5*fWidth*TMath::Tan(TMath::ACos(cosTheta1));
766 Double_t dL2 = 0.5*fWidth*TMath::Tan(TMath::ACos(cosTheta2));
767 if (localVect1[0]<0) dL1 = - dL1;
768 if (localVect2[0]<0) dL2 = - dL2;
770 Double_t cosPhi1 = -1./TMath::Sqrt( 1 + localVect1[1]*localVect1[1]
771 /localVect1[2]/localVect1[2] );
772 Double_t cosPhi2 = 1./TMath::Sqrt( 1 + localVect2[1]*localVect2[1]
773 /localVect2[2]/localVect2[2] );
774 if (localVect1[2]<0) cosPhi1 = -cosPhi1;
775 if (localVect2[2]<0) cosPhi2 = -cosPhi2;
777 Double_t tanACosCosPhi1 = TMath::Tan(TMath::ACos(cosPhi1));
778 Double_t tanACosCosPhi2 = TMath::Tan(TMath::ACos(cosPhi2));
779 if (localVect1[1]<0) tanACosCosPhi1 = -tanACosCosPhi1;
780 if (localVect2[1]<0) tanACosCosPhi2 = -tanACosCosPhi2;
782 Double_t dl1 = 0.5*fThick*tanACosCosPhi1*0.99999999999999;
783 Double_t dl2 = 0.5*fThick*tanACosCosPhi2*0.99999999999999;
784 // 0.9999999999999 is for correcting dawn problems in TGeo...
785 //=================================================
786 // Create the segment
787 TGeoArb8 *cableSeg = new TGeoArb8(fThick/2);
788 cableSeg->SetVertex( 0, -fWidth/2, -length/2 - dL1 + dl1);
789 cableSeg->SetVertex( 1, -fWidth/2, length/2 + dL2 - dl2);
790 cableSeg->SetVertex( 2, fWidth/2, length/2 - dL2 - dl2);
791 cableSeg->SetVertex( 3, fWidth/2, -length/2 + dL1 + dl1);
792 cableSeg->SetVertex( 4, -fWidth/2, -length/2 - dL1 - dl1);
793 cableSeg->SetVertex( 5, -fWidth/2, length/2 + dL2 + dl2);
794 cableSeg->SetVertex( 6, fWidth/2, length/2 - dL2 + dl2);
795 cableSeg->SetVertex( 7, fWidth/2, -length/2 + dL1 - dl1);
797 TGeoMedium *airSDD = gGeoManager->GetMedium("ITS_AIR$");
798 TGeoVolume *vCableSeg = new TGeoVolume(GetName(), cableSeg, airSDD);
800 // add all cable layers
801 for (Int_t iLay=0; iLay<fNlayer; iLay++) {
803 Double_t dl1Lay = 0.5*fLayThickness[iLay]*tanACosCosPhi1;
804 Double_t dl2Lay = 0.5*fLayThickness[iLay]*tanACosCosPhi2;
806 Double_t ztr = -fThick/2;
807 for (Int_t i=0;i<iLay; i++) ztr+= fLayThickness[i];
808 ztr+= fLayThickness[iLay]/2;
810 Double_t dl1LayS = ztr*tanACosCosPhi1;
811 Double_t dl2LayS = ztr*tanACosCosPhi2;
813 TGeoArb8 *lay = new TGeoArb8(fLayThickness[iLay]/2);
814 lay->SetVertex( 0, -fWidth/2, -length/2 - dL1 + dl1Lay - dl1LayS);
815 lay->SetVertex( 1, -fWidth/2, length/2 + dL2 - dl2Lay + dl2LayS);
816 lay->SetVertex( 2, fWidth/2, length/2 - dL2 - dl2Lay + dl2LayS);
817 lay->SetVertex( 3, fWidth/2, -length/2 + dL1 + dl1Lay - dl1LayS);
818 lay->SetVertex( 4, -fWidth/2, -length/2 - dL1 - dl1Lay - dl1LayS);
819 lay->SetVertex( 5, -fWidth/2, length/2 + dL2 + dl2Lay + dl2LayS);
820 lay->SetVertex( 6, fWidth/2, length/2 - dL2 + dl2Lay + dl2LayS);
821 lay->SetVertex( 7, fWidth/2, -length/2 + dL1 - dl1Lay - dl1LayS);
822 TGeoVolume *vLay = new TGeoVolume("vCableSegLay", lay, fLayMedia[iLay]);
823 vLay->SetLineColor(fLayColor[iLay]);
825 if (fTranslation[iLay]==0)
826 fTranslation[iLay] = new TGeoTranslation(0, 0, ztr);
827 vCableSeg->AddNode(vLay, iLay+1, fTranslation[iLay]);
830 vCableSeg->SetVisibility(kFALSE);
835 //________________________________________________________________________
836 TGeoVolume *AliITSv11GeomCableFlat::CreateCylSegment(Double_t &phi,
840 Double_t phi1 = 360-phi;
841 Double_t phi2 = 360+phi;
843 Double_t rMin = r-fThick/2;
844 Double_t rMax = r+fThick/2;
845 //=================================================
846 // Create the segment
848 TGeoTubeSeg *cableSeg = new TGeoTubeSeg(rMin, rMax, fWidth/2,
850 TGeoMedium *airSDD = gGeoManager->GetMedium("ITS_AIR$");
851 TGeoVolume *vCableSeg = new TGeoVolume(GetName(), cableSeg, airSDD);
853 // add all cable layers
854 for (Int_t iLay=0; iLay<fNlayer; iLay++) {
856 Double_t ztr = -fThick/2;
857 for (Int_t i=0;i<iLay; i++) ztr+= fLayThickness[i];
860 rMax = r + ztr + fLayThickness[iLay];
861 TGeoTubeSeg *lay = new TGeoTubeSeg(rMin, rMax, fWidth/2,
864 TGeoVolume *vLay = new TGeoVolume("vCableSegLay", lay, fLayMedia[iLay]);
865 vLay->SetLineColor(fLayColor[iLay]);
867 vCableSeg->AddNode(vLay, iLay+1, 0);
870 vCableSeg->SetVisibility(kFALSE);
875 //________________________________________________________________________
876 TGeoVolume *AliITSv11GeomCableFlat::CreateBoxSegment( Double_t *coord1,
880 //=================================================
881 // Create a segment for the case it is a simple box
882 Double_t dx = coord2[0]-coord1[0];
883 Double_t dy = coord2[1]-coord1[1];
884 Double_t dz = coord2[2]-coord1[2];
885 Double_t length = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
887 TGeoBBox *cableSeg = new TGeoBBox(fWidth/2, length/2, fThick/2);
889 TGeoMedium *airSDD = gGeoManager->GetMedium("ITS_AIR$");
890 TGeoVolume *vCableSeg = new TGeoVolume(GetName(), cableSeg, airSDD);
892 // add all cable layers
893 for (Int_t iLay=0; iLay<fNlayer; iLay++) {
895 Double_t ztr = -fThick/2;
896 for (Int_t i=0;i<iLay; i++) ztr+= fLayThickness[i];
897 ztr+= fLayThickness[iLay]/2;
899 TGeoBBox *lay = new TGeoBBox(fWidth/2, length/2, fLayThickness[iLay]/2);
902 TGeoVolume *vLay = new TGeoVolume("vCableSegLay", lay, fLayMedia[iLay]);
903 vLay->SetLineColor(fLayColor[iLay]);
905 if (fTranslation[iLay]==0)
906 fTranslation[iLay] = new TGeoTranslation(0, 0, ztr);
907 vCableSeg->AddNode(vLay, iLay+1, fTranslation[iLay]);
910 vCableSeg->SetVisibility(kFALSE);
914 //________________________________________________________________________
915 void AliITSv11GeomCableFlat::SetNLayers(Int_t nLayers) {
916 // Set the number of layers
917 if((nLayers>0) &&(nLayers<=fgkCableMaxLayer)) {
920 for (Int_t i=0; i<fgkCableMaxLayer ; i++) {
921 fLayThickness[i] = 0;
929 //________________________________________________________________________
930 Int_t AliITSv11GeomCableFlat::SetLayer(Int_t nLayer, Double_t thick,
931 TGeoMedium *medium, Int_t color) {
932 // Set the layer number nLayer
933 if ((nLayer<0)||(nLayer>=fNlayer)) {
934 printf("Set wrong layer number of the cable\n");
938 if (fLayThickness[nLayer-1]<=0) {
939 printf("AliITSv11GeomCableFlat::SetLayer():"
940 " You must define cable layer %i first !",nLayer-1);
944 Double_t thickTot = 0;
945 for (Int_t i=0; i<nLayer; i++) thickTot += fLayThickness[i];
947 if (thickTot-1e-10>fThick) {
948 printf("Can't add this layer, cable thickness would be higher than total\n");
952 fLayThickness[nLayer] = thick;
953 fLayMedia[nLayer] = medium;
954 fLayColor[nLayer] = color;
955 fTranslation[nLayer] = 0;