fix bugs add new features
[u/mrichter/AliRoot.git] / ITS / AliITSv11GeomCableRound.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
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  **************************************************************************/
15
16
17
18 // General Root includes
19 //#include <Riostream.h>
20 #include <TMath.h>
21 #include <TVectorD.h>
22
23 // Root Geometry includes
24 #include <TGeoManager.h>
25 #include <TGeoVolume.h>
26 #include <TGeoTube.h>
27 #include <TGeoTorus.h>
28 #include <TGeoMatrix.h>
29
30 #include "AliITSv11GeomCableRound.h"
31
32 //*************************************************************************
33 //   Class for round cables
34 //
35 // Ludovic Gaudichet                                   gaudichet@to.infn.it
36 //*************************************************************************
37
38 /*
39 // ************************************************************************
40 // Here is a example on how to use this class
41 // ************************************************************************
42
43   // Getting some media 
44   TGeoMedium *air   = gGeoManager->GetMedium("ITS_AIR$");
45   TGeoMedium *water = gGeoManager->GetMedium("ITS_WATER");
46   TGeoMedium *alu   = gGeoManager->GetMedium("ITS_ITSal"); 
47
48   // Creating a small box inside a bigger one (containers)
49   TGeoBBox *box1      = new TGeoBBox("box1", 6,10,10);
50   TGeoBBox *bigBox    = new TGeoBBox("bigBox", 20,10,10);
51   TGeoVolume *vbox1   = new TGeoVolume("vbox1", box1, air);
52   TGeoVolume *vBigBox = new TGeoVolume("vBigBox", bigBox, air);
53   vbox1->SetVisibility(kFALSE);
54   vBigBox->SetVisibility(kFALSE);
55
56   TGeoTranslation *tr1 = new TGeoTranslation("negTr",-14,0,0);
57   vBigBox->AddNode(vbox1, 1, tr1);
58   moth->AddNode(vBigBox, 1, 0);
59
60   // **************************************************
61   // Inserting a round cable (or here a water pipe...)
62   // **************************************************
63
64   Int_t waterColor = 7;
65   Int_t aluColor = 5;
66   AliITSv11GeomCableRound roundCable("waterPipe", 0.9); //radius of 0.9cm
67   roundCable.SetNLayers(2); 
68   roundCable.SetLayer(0, 0.7, water, waterColor); // radius of 0.7cm
69   roundCable.SetLayer(1, 0.2, alu, aluColor);     // thickness of 0.2cm
70
71   // ****** Set check points and their containers ******
72   // The 2 first points are in the small box (vbox1)
73   // The second point is at the boundary 
74
75   Double_t coord0[3] = {0,-2,-2};
76   Double_t coord1[3] = {6,2,1};
77   Double_t vect0[3]  = {1,1,0};
78   Double_t vect1[3]  = {1,0,0};
79   // coordinates have to be given in the specified container
80   // reference system (here it's going to be vbox1).
81   // vect1 and vect2 are vectors perpendicular to the segment ends
82   // (These vectors don't need to be normalized)
83   roundCable.AddCheckPoint( vbox1, 0, coord0, vect0);
84   roundCable.AddCheckPoint( vbox1, 1, coord1, vect1);
85
86   // Then, let's cross the boundary ! You just need
87   // to put the next point in the other volume, vBigBox.
88   // At the moment of creating the second segment, it will
89   // be inserted in this volume. That is why the point 1 had to
90   // be at the boundary, because otherwise the second segment
91   // between de points 1 and 2 would have been inserted in the
92   // vBigBox but in the same time would have cross its
93   // boundary ...
94   Double_t coord2[3] = {-2,6,4}; // coord. syst. of vBigBox !
95   Double_t vect2[3]= {1,1,0.5};
96   roundCable.AddCheckPoint( vBigBox, 2, coord2, vect2);
97
98   Double_t coord3[3] = {4,6,4};
99   Double_t vect3[3]= {-1,0,0};
100   roundCable.AddCheckPoint( vBigBox, 3, coord3, vect3);
101
102   Double_t coord4[3] = {4,0,-4};
103   Double_t vect4[3]= {1,0,0};
104   roundCable.AddCheckPoint( vBigBox, 4, coord4, vect4);
105
106   Double_t coord5[3] = {4,-6,4};
107   Double_t vect5[3]= {1,0,0};
108   roundCable.AddCheckPoint( vBigBox, 5, coord5, vect5);
109
110   Double_t coord6[3] = {7,-6,4};
111   Double_t vect6[3]= {1,0,0};
112   roundCable.AddCheckPoint( vBigBox, 6, coord6, vect6);
113
114   Double_t r = 7;
115   Double_t angle = 70*TMath::DegToRad(); 
116   Double_t coord7[3] = {coord6[0] +r*sin(angle), coord6[1],
117                         coord6[2] -r*(1-cos(angle)) };
118   Double_t vect7[3]= {r*cos(angle),0,-r*sin(angle)};
119   roundCable.AddCheckPoint( vBigBox, 7, coord7, vect7);
120
121   Double_t coord8[3] = { coord7[0]+vect7[0], coord7[1]+vect7[1],-10};
122   Double_t vect8[3]= {0,0,1};
123   roundCable.AddCheckPoint( vBigBox, 8, coord8, vect8);
124
125   // ****** Creating the corresponding volume ******
126   // Since the container volumes of the check points have
127   // been recorded, this can be done at any moments, providing
128   // that the container volumes are found in the sub-nodes
129   // of the initial node (the top volume of the TGeoManager or
130   // the volume set in SetInitialNode(TGeoVolume*) function)
131
132   roundCable.SetInitialNode(vBigBox); //Set the root node
133   roundCable.CreateAndInsertCableSegment( 1);
134   // This command means : create the segment between point 0
135   // and point 1. The segment is automatically inserted in the
136   // container volume of point 1.
137
138   roundCable.CreateAndInsertCableSegment( 2);
139   roundCable.CreateAndInsertCableSegment( 3);
140
141   // The following segment is going to be a torus segment.
142   // The radius and position of the torus is defined by the
143   // orthogonal vector of point 4 (the orientation of this vector
144   // and the position of the 2 check points are enough to define
145   // completely the torus)
146   roundCable.CreateAndInsertTorusSegment( 4, 180);
147   // The second argument is an additionnal rotation of the
148   // segment around the axis defined by the 2 check points.
149
150   roundCable.CreateAndInsertTorusSegment( 5);
151   roundCable.CreateAndInsertCableSegment( 6);
152   roundCable.CreateAndInsertTorusSegment( 7,180);
153   roundCable.CreateAndInsertCableSegment( 8);
154
155 */
156
157
158
159 ClassImp(AliITSv11GeomCableRound)
160
161 //________________________________________________________________________
162 AliITSv11GeomCableRound::
163 AliITSv11GeomCableRound(const char* name, Double_t radius) :
164   AliITSv11GeomCable(name),
165   fRadius(radius),
166   fNlayer(0),
167   fPhiMin(0),
168   fPhiMax(360)
169  {
170   // Constructor
171   for (Int_t i=0; i<fgkCableMaxLayer ; i++) {
172     fLayThickness[i] = 0;
173     fLayColor[i] = 0;
174     fLayMedia[i] = 0;
175   };
176 }
177 /*
178 //________________________________________________________________________
179 AliITSv11GeomCableRound::AliITSv11GeomCableRound(const AliITSv11GeomCableRound &s) :
180   AliITSv11GeomCable(s),fRadius(s.fRadius),fNlayer(s.fNlayer),fPhiMin(s.fPhiMin),
181   fPhiMax(s.fPhiMax)
182 {
183   //     Copy Constructor 
184   for (Int_t i=0; i<s.fNlayer; i++) {
185     fLayThickness[i] = s.fLayThickness[i];
186     fLayMedia[i] = s.fLayMedia[i];
187     fLayColor[i] = s.fLayColor[i];
188   }
189 }
190
191 //________________________________________________________________________
192 AliITSv11GeomCableRound& AliITSv11GeomCableRound::
193 operator=(const AliITSv11GeomCableRound &s) {
194   //     Assignment operator
195   // Not fully inplemented yet !!!
196
197   if(&s == this) return *this;
198   *this = s;
199   fRadius = s.fRadius;
200   fPhiMin = s.fPhiMin;
201   fPhiMax = s.fPhiMax;
202   fNlayer = s.fNlayer;
203   for (Int_t i=0; i<s.fNlayer; i++) {
204     fLayThickness[i] = s.fLayThickness[i];
205     fLayMedia[i] = s.fLayMedia[i];
206     fLayColor[i] = s.fLayColor[i];
207   };
208   return *this;
209 }
210 */
211 //________________________________________________________________________
212 Int_t AliITSv11GeomCableRound::GetPoint( Int_t iCheckPt, Double_t *coord)
213   const {
214   // Get check point #iCheckPt
215   TVectorD *coordVector =(TVectorD *)fPointArray.UncheckedAt(2*iCheckPt);
216 #if ROOT_VERSION_CODE < ROOT_VERSION(4,0,0)
217   CopyFrom(coord, coordVector->GetElements());
218 #else
219   CopyFrom(coord, coordVector->GetMatrixArray());
220 #endif 
221   return kTRUE;
222 }
223
224 //________________________________________________________________________
225 Int_t AliITSv11GeomCableRound::GetVect( Int_t iCheckPt, Double_t *coord)
226   const {
227   //
228   // Get vector transverse to the section at point #iCheckPt
229   //
230
231   TVectorD *coordVector =(TVectorD *)fPointArray.UncheckedAt(2*iCheckPt+1);
232 #if ROOT_VERSION_CODE < ROOT_VERSION(4,0,0)
233   CopyFrom(coord, coordVector->GetElements());
234 #else
235   CopyFrom(coord, coordVector->GetMatrixArray());
236 #endif 
237   return kTRUE;
238 }
239
240 //________________________________________________________________________
241 void AliITSv11GeomCableRound::AddCheckPoint( TGeoVolume *vol, Int_t iCheckPt,
242                                                Double_t *coord, Double_t *orthVect)
243 {
244   //
245   // Add point #iCheckPt and its transverse vector. Point is added at (i) in
246   // fPointArray and the vector is added at (i+1)
247   //
248
249
250   if (iCheckPt>=fVolumeArray.GetEntriesFast()) {
251     fVolumeArray.AddLast(vol);
252     TVectorD *point = new TVectorD(3,coord);
253     TVectorD *vect  = new TVectorD(3,orthVect);
254     fPointArray.AddLast(point);
255     fPointArray.AddLast(vect);
256
257   } else if ((iCheckPt >= 0)&&(iCheckPt < fVolumeArray.GetEntriesFast())) {
258     fVolumeArray.AddAt(vol, iCheckPt);
259     TVectorD *point = new TVectorD(3,coord);
260     TVectorD *vect  = new TVectorD(3,orthVect);
261     fPointArray.AddAt(point, iCheckPt*2  );
262     fPointArray.AddAt(vect,  iCheckPt*2+1);
263   };
264 }
265
266 //________________________________________________________________________
267 void AliITSv11GeomCableRound::PrintCheckPoints() const {
268   // Print all check points
269
270   printf("  ---\n  Printing all check points of the round cable\n");
271   for (Int_t i = 0; i<fVolumeArray.GetEntriesFast(); i++) {
272     TVectorD *coordVector = (TVectorD *)fPointArray.UncheckedAt(i*2);
273     //TVectorD *vectVector = (TVectorD *)fPointArray.UncheckedAt(i*2+1);
274     Double_t coord[3];
275 #if ROOT_VERSION_CODE < ROOT_VERSION(4,0,0)
276     CopyFrom(coord, coordVector->GetElements());
277 #else
278     CopyFrom(coord, coordVector->GetMatrixArray());
279 #endif 
280     printf("   ( %.2f, %.2f, %.2f )\n", coord[0], coord[1], coord[2]);
281   };
282 }
283
284
285 //________________________________________________________________________
286 TGeoVolume* AliITSv11GeomCableRound::CreateAndInsertCableSegment(Int_t p2,
287                                                                  TGeoCombiTrans** ct)
288 {
289 //    Creates a cable segment between points p1 and p2.
290 //
291 // The segment volume is created inside the volume containing point2
292 // Therefore this segment should be defined in this volume only.
293 // I mean here that, if the previous point is in another volume,
294 // it should be just at the border between the 2 volumes. Also the
295 // orientation vector of the previous point should be othogonal to
296 // the surface between the 2 volumes.
297
298   TGeoNode *mainNode;
299   if (fInitialNode==0) {
300     TObjArray *nodes = gGeoManager->GetListOfNodes();
301     if (nodes->GetEntriesFast()==0) return 0;
302     mainNode = (TGeoNode *) nodes->UncheckedAt(0);
303   } else {
304     mainNode = fInitialNode;
305   };
306
307   Int_t p1 = p2 - 1;
308   TGeoVolume *p1Vol = GetVolume(p1);
309   TGeoVolume *p2Vol = GetVolume(p2);
310
311   ResetCheckDaughter();
312   fCurrentVol = p1Vol;
313   if (! CheckDaughter(mainNode)) {
314     printf("Error::volume containing point is not visible in node tree!\n");
315     return 0;
316   };
317
318   Double_t coord1[3], coord2[3], vect1[3], vect2[3];
319   //=================================================
320   // Get p1 position in the systeme of p2
321   if (p1Vol!=p2Vol) {
322
323     Int_t p1nodeInd[fgkCableMaxNodeLevel]; 
324     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p1nodeInd[i]=fNodeInd[i];
325     Int_t p1volLevel = 0;
326     while (p1nodeInd[p1volLevel]!=-1) p1volLevel++;
327     p1volLevel--;
328
329     ResetCheckDaughter();
330     fCurrentVol = p2Vol;
331     if (! CheckDaughter(mainNode)) {
332       printf("Error::volume containing point is not visible in node tree!\n");
333       return 0;
334     };
335     Int_t p2nodeInd[fgkCableMaxNodeLevel];
336     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p2nodeInd[i]=fNodeInd[i];
337     Int_t commonMotherLevel = 0;
338     while (p1nodeInd[commonMotherLevel]==fNodeInd[commonMotherLevel])
339       commonMotherLevel++;
340     commonMotherLevel--;
341     Int_t p2volLevel = 0;
342     while (fNodeInd[p2volLevel]!=-1) p2volLevel++;
343     p2volLevel--;
344
345     // Get coord and vect of p1 in the common mother reference system
346     GetCheckPoint(p1, 0, p1volLevel-commonMotherLevel, coord1);
347     GetCheckVect( p1, 0, p1volLevel-commonMotherLevel, vect1);
348     // Translate them in the reference system of the volume containing p2    
349     TGeoNode *pathNode[fgkCableMaxNodeLevel];
350     pathNode[0] = mainNode;
351     for (Int_t i=0; i<=p2volLevel; i++) {
352       pathNode[i+1] = pathNode[i]->GetDaughter(p2nodeInd[i]);
353     };
354     Double_t globalCoord1[3] = {coord1[0], coord1[1], coord1[2]}; 
355     Double_t globalVect1[3]  = {vect1[0], vect1[1], vect1[2]};
356
357     for (Int_t i = commonMotherLevel+1; i<=p2volLevel; i++) {
358       pathNode[i+1]->GetMatrix()->MasterToLocal(globalCoord1, coord1);
359       pathNode[i+1]->GetMatrix()->MasterToLocalVect(globalVect1, vect1);
360       CopyFrom(globalCoord1, coord1);
361       CopyFrom(globalVect1, vect1);
362     };
363   } else {
364     GetCheckPoint(p1, 0, 0, coord1);
365     GetCheckVect(p1, 0, 0, vect1);
366   };
367   
368   //=================================================
369   // Get p2 position in the systeme of p2
370   GetCheckPoint(p2, 0, 0, coord2);
371   GetCheckVect(p2, 0, 0, vect2);
372
373   Double_t cx = (coord1[0]+coord2[0])/2;
374   Double_t cy = (coord1[1]+coord2[1])/2;
375   Double_t cz = (coord1[2]+coord2[2])/2;
376   Double_t dx = coord2[0]-coord1[0];
377   Double_t dy = coord2[1]-coord1[1];
378   Double_t dz = coord2[2]-coord1[2];
379
380   //=================================================
381   // Positionning of the segment between the 2 points
382   if ((dy<1e-31)&&(dy>0)) dy = 1e-31;
383   if ((dz<1e-31)&&(dz>0)) dz = 1e-31;
384   if ((dy>-1e-31)&&(dy<0)) dy = -1e-31;
385   if ((dz>-1e-31)&&(dz<0)) dz = -1e-31;
386
387   Double_t angleRot1 = -TMath::ATan2(dx,dy);
388   Double_t planDiagL = TMath::Sqrt(dy*dy+dx*dx);
389   Double_t angleRotDiag = -TMath::ATan2(planDiagL,dz);
390   TGeoRotation *rot = new TGeoRotation("",angleRot1*TMath::RadToDeg(),
391                                        angleRotDiag*TMath::RadToDeg(),
392                                        0);
393   Double_t localVect1[3], localVect2[3];
394   rot->MasterToLocalVect(vect1, localVect1);
395   rot->MasterToLocalVect(vect2, localVect2);
396   TGeoTranslation *trans = new TGeoTranslation("",cx, cy, cz);
397
398   //=================================================
399   // Create the segment and add it to the mother volume
400   TGeoVolume *vCableSeg = CreateSegment(coord1, coord2,
401                                         localVect1, localVect2, p2);
402
403   TGeoCombiTrans  *combi = new TGeoCombiTrans(*trans, *rot);
404   p2Vol->AddNode(vCableSeg, p2, combi);
405   //=================================================
406   delete rot;
407   delete trans;
408
409   if (fDebug) {
410     printf("---\n  Cable segment points : ");
411     printf("%f, %f, %f\n",coord1[0], coord1[1], coord1[2]);
412     printf("%f, %f, %f\n",coord2[0], coord2[1], coord2[2]);
413   };
414 //   #include <TGeoSphere.h>
415 //   TGeoMedium *airSDD = gGeoManager->GetMedium("ITS_AIR$");
416 //   TGeoSphere *sphere = new TGeoSphere(0, 0.15);
417 //   TGeoVolume *vSphere = new TGeoVolume("", sphere, airSDD);
418 //   TGeoTranslation *trC = new TGeoTranslation("", cx, cy, cz);
419 //   TGeoTranslation *tr1 = new TGeoTranslation("",coord1[0],
420 //                                           coord1[1],coord1[2]);
421 //   TGeoTranslation *tr2 = new TGeoTranslation("",coord2[0],
422 //                                           coord2[1],coord2[2]);
423 //   p2Vol->AddNode(vSphere, p2*3-2, trC);
424 //   p2Vol->AddNode(vSphere, p2*3-1, tr1);
425 //   p2Vol->AddNode(vSphere, p2*3  , tr2);
426
427   if (ct) *ct = combi;
428   return vCableSeg;
429 }
430
431 //________________________________________________________________________
432 TGeoVolume* AliITSv11GeomCableRound::CreateAndInsertTubeSegment(Int_t p2,
433                                                                  TGeoCombiTrans** ct)
434 {
435 //    Creates a cable segment between points p1 and p2.
436 //
437 //  This creates simple tube sections, i.e. the cable ends are
438 // cutted perpendicularly to the tube axis. The method has to
439 // be used only in this simple case, in ordder to save some memory
440
441   TGeoNode *mainNode;
442   if (fInitialNode==0) {
443     TObjArray *nodes = gGeoManager->GetListOfNodes();
444     if (nodes->GetEntriesFast()==0) return 0;
445     mainNode = (TGeoNode *) nodes->UncheckedAt(0);
446   } else {
447     mainNode = fInitialNode;
448   };
449
450   Int_t p1 = p2 - 1;
451   TGeoVolume *p1Vol = GetVolume(p1);
452   TGeoVolume *p2Vol = GetVolume(p2);
453
454   ResetCheckDaughter();
455   fCurrentVol = p1Vol;
456   if (! CheckDaughter(mainNode)) {
457     printf("Error::volume containing point is not visible in node tree!\n");
458     return 0;
459   };
460
461   Double_t coord1[3], coord2[3], vect1[3], vect2[3];
462   //=================================================
463   // Get p1 position in the systeme of p2
464   if (p1Vol!=p2Vol) {
465
466     Int_t p1nodeInd[fgkCableMaxNodeLevel]; 
467     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p1nodeInd[i]=fNodeInd[i];
468     Int_t p1volLevel = 0;
469     while (p1nodeInd[p1volLevel]!=-1) p1volLevel++;
470     p1volLevel--;
471
472     ResetCheckDaughter();
473     fCurrentVol = p2Vol;
474     if (! CheckDaughter(mainNode)) {
475       printf("Error::volume containing point is not visible in node tree!\n");
476       return 0;
477     };
478     Int_t p2nodeInd[fgkCableMaxNodeLevel];
479     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p2nodeInd[i]=fNodeInd[i];
480     Int_t commonMotherLevel = 0;
481     while (p1nodeInd[commonMotherLevel]==fNodeInd[commonMotherLevel])
482       commonMotherLevel++;
483     commonMotherLevel--;
484     Int_t p2volLevel = 0;
485     while (fNodeInd[p2volLevel]!=-1) p2volLevel++;
486     p2volLevel--;
487
488     // Get coord and vect of p1 in the common mother reference system
489     GetCheckPoint(p1, 0, p1volLevel-commonMotherLevel, coord1);
490     GetCheckVect( p1, 0, p1volLevel-commonMotherLevel, vect1);
491     // Translate them in the reference system of the volume containing p2    
492     TGeoNode *pathNode[fgkCableMaxNodeLevel];
493     pathNode[0] = mainNode;
494     for (Int_t i=0; i<=p2volLevel; i++) {
495       pathNode[i+1] = pathNode[i]->GetDaughter(p2nodeInd[i]);
496     };
497     Double_t globalCoord1[3] = {coord1[0], coord1[1], coord1[2]}; 
498     Double_t globalVect1[3]  = {vect1[0], vect1[1], vect1[2]};
499
500     for (Int_t i = commonMotherLevel+1; i<=p2volLevel; i++) {
501       pathNode[i+1]->GetMatrix()->MasterToLocal(globalCoord1, coord1);
502       pathNode[i+1]->GetMatrix()->MasterToLocalVect(globalVect1, vect1);
503       CopyFrom(globalCoord1, coord1);
504       CopyFrom(globalVect1, vect1);
505     };
506   } else {
507     GetCheckPoint(p1, 0, 0, coord1);
508     GetCheckVect(p1, 0, 0, vect1);
509   };
510   
511   //=================================================
512   // Get p2 position in the systeme of p2
513   GetCheckPoint(p2, 0, 0, coord2);
514   GetCheckVect(p2, 0, 0, vect2);
515
516   Double_t cx = (coord1[0]+coord2[0])/2;
517   Double_t cy = (coord1[1]+coord2[1])/2;
518   Double_t cz = (coord1[2]+coord2[2])/2;
519   Double_t dx = coord2[0]-coord1[0];
520   Double_t dy = coord2[1]-coord1[1];
521   Double_t dz = coord2[2]-coord1[2];
522
523   //=================================================
524   // Positionning of the segment between the 2 points
525   if ((dy<1e-31)&&(dy>0)) dy = 1e-31;
526   if ((dz<1e-31)&&(dz>0)) dz = 1e-31;
527   if ((dy>-1e-31)&&(dy<0)) dy = -1e-31;
528   if ((dz>-1e-31)&&(dz<0)) dz = -1e-31;
529
530   Double_t angleRot1 = -TMath::ATan2(dx,dy);
531   Double_t planDiagL = TMath::Sqrt(dy*dy+dx*dx);
532   Double_t angleRotDiag = -TMath::ATan2(planDiagL,dz);
533   TGeoRotation *rot = new TGeoRotation("",angleRot1*TMath::RadToDeg(),
534                                        angleRotDiag*TMath::RadToDeg(),
535                                        0);
536   TGeoTranslation *trans = new TGeoTranslation("",cx, cy, cz);
537
538   //=================================================
539   // Create the segment and add it to the mother volume
540   TGeoVolume *vCableSeg = CreateTubeSegment( coord1,coord2, p2);
541
542   TGeoCombiTrans  *combi = new TGeoCombiTrans(*trans, *rot);
543   p2Vol->AddNode(vCableSeg, p2, combi);
544   //=================================================
545   delete rot;
546   delete trans;
547
548   if (fDebug) {
549     printf("---\n  Cable segment points : ");
550     printf("%f, %f, %f\n",coord1[0], coord1[1], coord1[2]);
551     printf("%f, %f, %f\n",coord2[0], coord2[1], coord2[2]);
552   };
553
554   if (ct) *ct = combi;
555   return vCableSeg;
556 }
557
558 //________________________________________________________________________
559 TGeoVolume* AliITSv11GeomCableRound::CreateAndInsertTorusSegment(Int_t p2,
560                                                                  Double_t rotation,
561                                                                  TGeoCombiTrans** ct)
562 {
563   // Create a torus cable segment between points p1 and p2.
564   // The radius and position of the torus is defined by the
565   // perpendicular vector of point p2 (the orientation of this vector
566   // and the position of the 2 check points are enough to completely
567   // define the torus)
568
569   TGeoNode *mainNode;
570   if (fInitialNode==0) {
571     TObjArray *nodes = gGeoManager->GetListOfNodes();
572     if (nodes->GetEntriesFast()==0) return 0;
573     mainNode = (TGeoNode *) nodes->UncheckedAt(0);
574   } else {
575     mainNode = fInitialNode;
576   };
577
578   Int_t p1 = p2 - 1;
579   TGeoVolume *p1Vol = GetVolume(p1);
580   TGeoVolume *p2Vol = GetVolume(p2);
581
582   ResetCheckDaughter();
583   fCurrentVol = p1Vol;
584   if (! CheckDaughter(mainNode)) {
585     printf("Error::volume containing point is not visible in node tree!\n");
586     return 0;
587   };
588
589   Double_t coord1[3], coord2[3], vect1[3], vect2[3];
590   //=================================================
591   // Get p1 position in the systeme of p2
592   if (p1Vol!=p2Vol) {
593
594     Int_t p1nodeInd[fgkCableMaxNodeLevel]; 
595     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p1nodeInd[i]=fNodeInd[i];
596     Int_t p1volLevel = 0;
597     while (p1nodeInd[p1volLevel]!=-1) p1volLevel++;
598     p1volLevel--;
599
600     ResetCheckDaughter();
601     fCurrentVol = p2Vol;
602     if (! CheckDaughter(mainNode)) {
603       printf("Error::volume containing point is not visible in node tree!\n");
604       return 0;
605     };
606     Int_t p2nodeInd[fgkCableMaxNodeLevel];
607     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p2nodeInd[i]=fNodeInd[i];
608     Int_t commonMotherLevel = 0;
609     while (p1nodeInd[commonMotherLevel]==fNodeInd[commonMotherLevel])
610       commonMotherLevel++;
611     commonMotherLevel--;
612     Int_t p2volLevel = 0;
613     while (fNodeInd[p2volLevel]!=-1) p2volLevel++;
614     p2volLevel--;
615
616     // Get coord and vect of p1 in the common mother reference system
617     GetCheckPoint(p1, 0, p1volLevel-commonMotherLevel, coord1);
618     GetCheckVect( p1, 0, p1volLevel-commonMotherLevel, vect1);
619     // Translate them in the reference system of the volume containing p2    
620     TGeoNode *pathNode[fgkCableMaxNodeLevel];
621     pathNode[0] = mainNode;
622     for (Int_t i=0; i<=p2volLevel; i++) {
623       pathNode[i+1] = pathNode[i]->GetDaughter(p2nodeInd[i]);
624     };
625     Double_t globalCoord1[3] = {coord1[0], coord1[1], coord1[2]}; 
626     Double_t globalVect1[3]  = {vect1[0], vect1[1], vect1[2]};
627
628     for (Int_t i = commonMotherLevel+1; i<=p2volLevel; i++) {
629       pathNode[i+1]->GetMatrix()->MasterToLocal(globalCoord1, coord1);
630       pathNode[i+1]->GetMatrix()->MasterToLocalVect(globalVect1, vect1);
631       CopyFrom(globalCoord1, coord1);
632       CopyFrom(globalVect1, vect1);
633     };
634   } else {
635     GetCheckPoint(p1, 0, 0, coord1);
636     GetCheckVect(p1, 0, 0, vect1);
637   };
638   
639   //=================================================
640   // Get p2 position in the systeme of p2
641   GetCheckPoint(p2, 0, 0, coord2);
642   GetCheckVect(p2, 0, 0, vect2);
643
644   Double_t cx = (coord1[0]+coord2[0])/2;
645   Double_t cy = (coord1[1]+coord2[1])/2;
646   Double_t cz = (coord1[2]+coord2[2])/2;
647   Double_t dx = coord2[0]-coord1[0];
648   Double_t dy = coord2[1]-coord1[1];
649   Double_t dz = coord2[2]-coord1[2];
650   Double_t length = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
651
652   //=================================================
653   // Positionning of the segment between the 2 points
654   if ((dy<1e-31)&&(dy>0)) dy = 1e-31;
655   if ((dz<1e-31)&&(dz>0)) dz = 1e-31;
656   if ((dy>-1e-31)&&(dy<0)) dy = -1e-31;
657   if ((dz>-1e-31)&&(dz<0)) dz = -1e-31;
658
659   Double_t angleRot1 = -TMath::ATan2(dx,dy);
660   Double_t planDiagL = TMath::Sqrt(dy*dy+dx*dx);
661   Double_t angleRotDiag = -TMath::ATan2(planDiagL,dz);
662
663   TGeoRotation rotTorusTemp("",angleRot1*TMath::RadToDeg(),
664                             angleRotDiag*TMath::RadToDeg(),0);
665   TGeoRotation rotTorusToZ("",0,90,0);
666   rotTorusTemp.MultiplyBy(&rotTorusToZ, kTRUE);
667   Double_t localVect2[3];
668   rotTorusTemp.MasterToLocalVect(vect2, localVect2);
669   if (localVect2[1]<0) {
670     localVect2[0] = -localVect2[0];
671     localVect2[1] = -localVect2[1];
672     localVect2[2] = -localVect2[2];
673   };
674   Double_t normVect2 = TMath::Sqrt(localVect2[0]*localVect2[0]+
675                                    localVect2[1]*localVect2[1]+
676                                    localVect2[2]*localVect2[2]);
677   Double_t axisX[3] = {1,0,0};
678   Double_t cosangleTorusSeg = (localVect2[0]*axisX[0]+
679                                localVect2[1]*axisX[1]+
680                                localVect2[2]*axisX[2])/normVect2;
681   Double_t angleTorusSeg = TMath::ACos(cosangleTorusSeg)*TMath::RadToDeg();
682
683   TGeoRotation rotTorus("",angleRot1*TMath::RadToDeg(),
684                         angleRotDiag*TMath::RadToDeg(),
685                         180-angleTorusSeg+rotation);
686   rotTorus.MultiplyBy(&rotTorusToZ, kTRUE);
687   rotTorus.MasterToLocalVect(vect2, localVect2);
688   if (localVect2[1]<0) {
689     localVect2[0] = -localVect2[0];
690     localVect2[1] = -localVect2[1];
691     localVect2[2] = -localVect2[2];
692   };
693   normVect2 = TMath::Sqrt(localVect2[0]*localVect2[0]+
694                           localVect2[1]*localVect2[1]+
695                           localVect2[2]*localVect2[2]);
696   Double_t axisY[3] = {0,1,0};
697   Double_t cosPhi = (localVect2[0]*axisY[0]+localVect2[1]*axisY[1]+
698                      localVect2[2]*axisY[2])/normVect2;
699   Double_t torusPhi1 = TMath::ACos(cosPhi);
700   Double_t torusR = (length/2)/TMath::Sin(torusPhi1);
701   torusPhi1 = torusPhi1*TMath::RadToDeg();
702   Double_t perpLength = TMath::Sqrt(torusR*torusR-length*length/4);
703   Double_t localTransT[3] = {-perpLength,0,0};
704   Double_t globalTransT[3];
705   rotTorus.LocalToMasterVect(localTransT, globalTransT);
706   TGeoTranslation transTorus("",cx+globalTransT[0],cy+globalTransT[1],
707                              cz+globalTransT[2]);
708
709   TGeoCombiTrans  *combiTorus = new TGeoCombiTrans(transTorus, rotTorus);
710
711   //=================================================
712   // Create the segment and add it to the mother volume
713   TGeoVolume *vCableSegT = CreateTorus(torusPhi1, torusR, p2);
714   p2Vol->AddNode(vCableSegT, p2, combiTorus);
715
716   if (fDebug) {
717     printf("---\n  Cable segment points : ");
718     printf("%f, %f, %f\n",coord1[0], coord1[1], coord1[2]);
719     printf("%f, %f, %f\n",coord2[0], coord2[1], coord2[2]);
720   };
721
722   if (ct) *ct = combiTorus;
723   return vCableSegT;
724 }
725
726 //________________________________________________________________________
727 TGeoVolume *AliITSv11GeomCableRound::CreateSegment( const Double_t *coord1,
728                                                     const Double_t *coord2,
729                                                       Double_t *localVect1,
730                                                       Double_t *localVect2, Int_t p)
731 {
732   // Create a cylindrical segment and its layers. The tube section is cutted by
733   // two planes, defined by the normal vectors localVect1 and localVect2
734
735   //=================================================
736   // Calculate segment "deformation"
737   Double_t dx = coord2[0]-coord1[0];
738   Double_t dy = coord2[1]-coord1[1];
739   Double_t dz = coord2[2]-coord1[2];
740   Double_t length = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
741
742   // normal vectors have to point outside the TGeoCtub :
743   if (-localVect1[2]<0) {
744     localVect1[0] = -localVect1[0];
745     localVect1[1] = -localVect1[1];
746     localVect1[2] = -localVect1[2];
747   };
748   if (localVect2[2]<0) {
749     localVect2[0] = -localVect2[0];
750     localVect2[1] = -localVect2[1];
751     localVect2[2] = -localVect2[2];
752   };
753   //=================================================
754   // Create the segment
755   TGeoCtub *cableSeg = new TGeoCtub(0, fRadius, length/2, fPhiMin, fPhiMax,
756                                     localVect1[0],localVect1[1],localVect1[2],
757                                     localVect2[0],localVect2[1],localVect2[2]);
758
759   TGeoMedium *skinMedia = fLayMedia[fNlayer-1];
760   char name[100];
761   snprintf(name, 100, "%s_%i",GetName(), p);
762   TGeoVolume *vCableSeg = new TGeoVolume(name, cableSeg, skinMedia);
763   vCableSeg->SetLineColor(fLayColor[fNlayer-1]);
764
765   // add all cable layers
766   Double_t layThickness[100+1];                        // 100 layers max !!!
767   layThickness[0] = 0;
768   for (Int_t iLay=0; iLay<fNlayer-1; iLay++) {
769     
770     layThickness[iLay+1] = fLayThickness[iLay]+layThickness[iLay];
771     TGeoCtub *lay = new TGeoCtub(layThickness[iLay], layThickness[iLay+1],
772                                  length/2, fPhiMin, fPhiMax,
773                                  localVect1[0],localVect1[1],localVect1[2],
774                                  localVect2[0],localVect2[1],localVect2[2]);
775
776     TGeoVolume *vLay = new TGeoVolume("vCableSegLay", lay, fLayMedia[iLay]);
777     vLay->SetLineColor(fLayColor[iLay]);
778     vCableSeg->AddNode(vLay, iLay+1, 0);
779   };
780
781   //vCableSeg->SetVisibility(kFALSE);
782   return vCableSeg;
783 }
784
785
786 //________________________________________________________________________
787 TGeoVolume *AliITSv11GeomCableRound::CreateTubeSegment( const Double_t *coord1,
788                                                         const Double_t *coord2,
789                                                         Int_t p)
790 {
791   // Create a cylindrical segment and its layers
792
793   //=================================================
794   // Calculate segment "deformation"
795   Double_t dx = coord2[0]-coord1[0];
796   Double_t dy = coord2[1]-coord1[1];
797   Double_t dz = coord2[2]-coord1[2];
798   Double_t length = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
799
800  //=================================================
801   // Create the segment
802
803   TGeoTubeSeg *cableSeg = new TGeoTubeSeg(0, fRadius, length/2, fPhiMin, fPhiMax);
804
805   TGeoMedium *skinMedia = fLayMedia[fNlayer-1];
806   char name[100];
807   snprintf(name, 100, "%s_%i",GetName(), p);
808   TGeoVolume *vCableSeg = new TGeoVolume(name, cableSeg, skinMedia);
809   vCableSeg->SetLineColor(fLayColor[fNlayer-1]);
810
811   // add all cable layers
812   Double_t layThickness[100+1];                        // 100 layers max !!!
813   layThickness[0] = 0;
814   for (Int_t iLay=0; iLay<fNlayer-1; iLay++) {
815     
816     layThickness[iLay+1] = fLayThickness[iLay]+layThickness[iLay];
817     TGeoTubeSeg*lay = new TGeoTubeSeg(layThickness[iLay], layThickness[iLay+1],
818                                       length/2, fPhiMin, fPhiMax);
819     TGeoVolume *vLay = new TGeoVolume("vCableSegLay", lay, fLayMedia[iLay]);
820     vLay->SetLineColor(fLayColor[iLay]);
821     vCableSeg->AddNode(vLay, iLay+1, 0);
822   };
823
824   //vCableSeg->SetVisibility(kFALSE);
825   return vCableSeg;
826 }
827
828
829 //________________________________________________________________________
830 TGeoVolume *AliITSv11GeomCableRound::CreateTorus( const Double_t &phi,
831                                                   const Double_t &r, Int_t p)
832 {
833   // Create one torus segment and its layers
834
835   Double_t torusR = r;
836 //   Double_t torusPhi1 = phi;
837 //   Double_t torusDPhi = -2*torusPhi1;  // bug in root ...
838   Double_t torusPhi1 = 360-phi;
839   Double_t torusDPhi = 2*phi;
840
841   //  // Create the segment, it will also work as the last layer
842   TGeoTorus *cableSeg = new TGeoTorus(torusR, 0, fRadius, torusPhi1, torusDPhi);
843   TGeoMedium *skinMedia = fLayMedia[fNlayer-1];
844   char name[100];
845   snprintf(name, 100, "%s_%i",GetName(),p);
846   TGeoVolume *vCableSeg = new TGeoVolume(name, cableSeg, skinMedia);
847   vCableSeg->SetLineColor(fLayColor[fNlayer-1]);
848
849   // add all cable layers but last
850   Double_t layThickness[100+1];                        // 100 layers max !!!
851   layThickness[0] = 0;
852   for (Int_t iLay=0; iLay<fNlayer-1; iLay++) {
853     
854     layThickness[iLay+1] = fLayThickness[iLay]+layThickness[iLay];
855     TGeoTorus *lay = new TGeoTorus(torusR, layThickness[iLay],
856                                    layThickness[iLay+1],
857                                    torusPhi1,torusDPhi);
858
859     TGeoVolume *vLay = new TGeoVolume("vCableSegLay",lay,fLayMedia[iLay]);
860     vLay->SetLineColor(fLayColor[iLay]);
861
862     vCableSeg->AddNode(vLay, iLay+1,0);
863   };
864
865   //vCableSeg->SetVisibility(kFALSE);
866   return vCableSeg;
867 }
868
869 //________________________________________________________________________
870 void AliITSv11GeomCableRound::SetNLayers(Int_t nLayers) {
871   // Set the total number of layers
872   if((nLayers>0) &&(nLayers<=fgkCableMaxLayer)) {
873     fNlayer = nLayers;
874     for (Int_t i = 0; i<fNlayer; i++) {
875       fLayThickness[i] = 0;
876       fLayMedia[i] = 0;
877     };
878   };
879 }
880
881 //________________________________________________________________________
882 Int_t AliITSv11GeomCableRound::SetLayer(Int_t nLayer, Double_t thick,
883                                            TGeoMedium *medium, Int_t color) {
884   // Set layer #nLayer
885   if ((nLayer<0)||(nLayer>=fNlayer)) {
886     printf("Set wrong layer number of the cable\n");
887     return kFALSE;
888   };
889   if (nLayer>0)
890     if (fLayThickness[nLayer-1]<=0) {
891       printf("You must define cable layer %i first !",nLayer-1);
892       return kFALSE;
893     };
894
895   Double_t thickTot = 0;
896   for (Int_t i=0; i<nLayer; i++) thickTot += fLayThickness[i];
897   thickTot += thick;
898   if (thickTot-1e-10>fRadius) {
899     printf("Can't add this layer, cable thickness would be higher than total\n");
900     return kFALSE;
901   };
902
903   fLayThickness[nLayer] = thick;
904   fLayMedia[nLayer] = medium;
905   fLayColor[nLayer] = color;
906
907   return kTRUE;
908 }