Changing once more (hopefully we get it correct this time...) the logic to trig the...
[u/mrichter/AliRoot.git] / ITS / AliITSv11GeomCableFlat.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 // $Id$
17
18 //*************************************************************************
19 //   Class for flat cables
20 //
21 // Ludovic Gaudichet                                   gaudichet@to.infn.it
22 //*************************************************************************
23
24
25
26 // General Root includes
27 //#include <Riostream.h>
28 #include <TMath.h>
29 #include <TVectorD.h>
30
31 // Root Geometry includes
32 #include <TGeoManager.h>
33 #include <TGeoVolume.h>
34 #include <TGeoArb8.h>
35 #include <TGeoTube.h>
36 #include <TGeoMatrix.h>
37 #include <TGeoNode.h>
38
39 #include "AliITSv11GeomCableFlat.h"
40
41
42 ClassImp(AliITSv11GeomCableFlat)
43
44 //________________________________________________________________________
45 AliITSv11GeomCableFlat::AliITSv11GeomCableFlat():
46   AliITSv11GeomCable(),
47   fWidth(0),
48   fThick(0),
49   fNlayer(0)
50 {
51   // constructor
52   for (Int_t i=0; i<fgkCableMaxLayer ; i++) {
53     fLayThickness[i] = 0;
54     fTranslation[i]  = 0;
55     fLayColor[i]     = 0;
56     fLayMedia[i]     = 0;  
57  };
58 }
59
60 //________________________________________________________________________
61 AliITSv11GeomCableFlat::
62 AliITSv11GeomCableFlat(const char* name, Double_t width, Double_t thick) :
63   AliITSv11GeomCable(name),
64   fWidth(width),
65   fThick(thick),
66   fNlayer(0)
67  {
68   // standard constructor
69   for (Int_t i=0; i<fgkCableMaxLayer ; i++) {
70     fLayThickness[i] = 0;
71     fTranslation[i]  = 0;
72     fLayColor[i]     = 0;
73     fLayMedia[i]     = 0;  
74   }; 
75 }
76
77 //________________________________________________________________________
78 AliITSv11GeomCableFlat::AliITSv11GeomCableFlat(const AliITSv11GeomCableFlat &s) :
79   AliITSv11GeomCable(s),fWidth(s.fWidth),fThick(s.fThick),fNlayer(s.fNlayer)
80 {
81   //     Copy Constructor 
82   for (Int_t i=0; i<s.fNlayer; i++) {
83     fLayThickness[i] = s.fLayThickness[i];
84     fTranslation[i] = s.fTranslation[i];
85     fLayMedia[i] = s.fLayMedia[i];
86     fLayColor[i] = s.fLayColor[i];
87   }
88 }
89
90 //________________________________________________________________________
91 AliITSv11GeomCableFlat& AliITSv11GeomCableFlat::
92 operator=(const AliITSv11GeomCableFlat &s) {
93   //     Assignment operator
94   // Not fully inplemented yet !!!
95
96   if(&s == this) return *this;
97   *this = s;
98   fWidth = s.fWidth;
99   fThick = s.fThick;
100   fNlayer = s.fNlayer;
101   for (Int_t i=0; i<s.fNlayer; i++) {
102     fLayThickness[i] = s.fLayThickness[i];
103     fTranslation[i] = s.fTranslation[i];
104     fLayMedia[i] = s.fLayMedia[i];
105     fLayColor[i] = s.fLayColor[i];
106   };
107   return *this;
108 }
109
110 //________________________________________________________________________
111 Int_t AliITSv11GeomCableFlat::GetPoint( Int_t iCheckPt, Double_t *coord)
112   const {
113   // Get the correct point #iCheckPt
114   TVectorD *coordVector =(TVectorD *)fPointArray.At(2*iCheckPt);
115   if (coordVector) {
116 #if ROOT_VERSION_CODE < ROOT_VERSION(4,0,0)
117     CopyFrom(coord, coordVector->GetElements());
118 #else
119     CopyFrom(coord, coordVector->GetMatrixArray());
120 #endif 
121     return kTRUE;
122   } else {
123     return kFALSE;
124   };
125 }
126
127 //________________________________________________________________________
128 Int_t AliITSv11GeomCableFlat::GetVect( Int_t iCheckPt, Double_t *coord)
129   const {
130   // Get the correct vect corresponding to point #iCheckPt
131
132   TVectorD *coordVector =(TVectorD *)fPointArray.At(2*iCheckPt+1);
133   if (coordVector) {
134 #if ROOT_VERSION_CODE < ROOT_VERSION(4,0,0)
135     CopyFrom(coord, coordVector->GetElements());
136 #else
137     CopyFrom(coord, coordVector->GetMatrixArray());
138 #endif 
139     return kTRUE;
140   } else {
141     return kFALSE;
142   };
143 }
144
145 //________________________________________________________________________
146 void AliITSv11GeomCableFlat::AddCheckPoint( TGeoVolume *vol, Int_t iCheckPt,
147                                             Double_t *coord, Double_t *orthVect)
148 {
149   //
150   // Add a check point. In the fPointArray, the point is at i and its vector
151   // is at i+1.
152   //
153
154 //   if (iCheckPt>=fVolumeArray.GetEntriesFast()) {
155 //     fVolumeArray.AddLast(vol);
156 //     TVectorD *point = new TVectorD(3,coord);
157 //     TVectorD *vect  = new TVectorD(3,orthVect);
158 //     fPointArray.AddLast(point);
159 //     fPointArray.AddLast(vect);
160
161 //   } else if ((iCheckPt >= 0)&&(iCheckPt < fVolumeArray.GetEntriesFast())) {
162 //     fVolumeArray.AddAt(vol, iCheckPt);
163 //     TVectorD *point = new TVectorD(3,coord);
164 //     TVectorD *vect  = new TVectorD(3,orthVect);
165 //     fPointArray.AddAt(point, iCheckPt*2  );
166 //     fPointArray.AddAt(vect,  iCheckPt*2+1);
167 //   };
168   fVolumeArray.AddAtAndExpand(vol, iCheckPt);
169   TVectorD *point = new TVectorD(3,coord);
170   TVectorD *vect  = new TVectorD(3,orthVect);
171   fPointArray.AddAtAndExpand(point, iCheckPt*2  );
172   fPointArray.AddAtAndExpand(vect,  iCheckPt*2+1);
173 }
174
175 //________________________________________________________________________
176 void AliITSv11GeomCableFlat::PrintCheckPoints() const {
177   // print all check points of the cable
178   printf("  ---\n  Printing all check points of the flat cable\n");
179   for (Int_t i = 0; i<fVolumeArray.GetEntriesFast(); i++) {
180      Double_t coord[3];
181      if (GetPoint( i, coord))
182        printf("   ( %.2f, %.2f, %.2f )\n", coord[0], coord[1], coord[2]);
183   };
184 }
185
186 //________________________________________________________________________
187 TGeoVolume* AliITSv11GeomCableFlat::CreateAndInsertCableSegment(Int_t p2,
188                                                                 Double_t rotation,
189                                                                 TGeoCombiTrans** ct)
190 {
191 //    Creates a cable segment between points p1 and p2.
192 //    Rotation is the eventual rotation of the flat cable
193 //    along its length axis
194 //
195 // The segment volume is created inside the volume containing point2
196 // Therefore this segment should be defined in this volume only.
197 // I mean here that, if the previous point is in another volume,
198 // it should be just at the border between the 2 volumes. Also the
199 // orientation vector of the previous point should be orthogonal to
200 // the surface between the 2 volumes.
201
202   TGeoNode *mainNode;
203   if (fInitialNode==0) {
204     TObjArray *nodes = gGeoManager->GetListOfNodes();
205     if (nodes->GetEntriesFast()==0) return 0;
206     mainNode = (TGeoNode *) nodes->UncheckedAt(0);
207   } else {
208     mainNode = fInitialNode;
209   };
210
211   Int_t p1 = p2 - 1;
212   TGeoVolume *p2Vol = GetVolume(p2);
213   TGeoVolume *p1Vol = GetVolume(p1);
214
215   ResetCheckDaughter();
216   fCurrentVol = p1Vol;
217   if (! CheckDaughter(mainNode)) {
218     printf("Error::volume containing point is not visible in node tree!\n");
219     return 0;
220   };
221
222   Double_t coord1[3], coord2[3], vect1[3], vect2[3];
223   //=================================================
224   // Get p1 position in the systeme of p2
225   if (p1Vol!=p2Vol) {
226
227     Int_t p1nodeInd[fgkCableMaxNodeLevel]; 
228     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p1nodeInd[i]=fNodeInd[i];
229     Int_t p1volLevel = 0;
230     while (p1nodeInd[p1volLevel]!=-1) p1volLevel++;
231     p1volLevel--;
232
233     ResetCheckDaughter();
234     fCurrentVol = p2Vol;
235     if (! CheckDaughter(mainNode)) {
236       printf("Error::volume containing point is not visible in node tree!\n");
237       return 0;
238     };
239     Int_t p2nodeInd[fgkCableMaxNodeLevel];
240     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p2nodeInd[i]=fNodeInd[i];
241     Int_t commonMotherLevel = 0;
242     while (p1nodeInd[commonMotherLevel]==fNodeInd[commonMotherLevel])
243       commonMotherLevel++;
244     commonMotherLevel--;
245     Int_t p2volLevel = 0;
246     while (fNodeInd[p2volLevel]!=-1) p2volLevel++;
247     p2volLevel--;
248
249     // Get coord and vect of p1 in the common mother reference system
250     if (! GetCheckPoint(p1, 0, p1volLevel-commonMotherLevel, coord1) )
251       return 0;
252     if (! GetCheckVect( p1, 0, p1volLevel-commonMotherLevel, vect1) )
253       return 0;
254
255     // Translate them in the reference system of the volume containing p2    
256     TGeoNode *pathNode[fgkCableMaxNodeLevel];
257     pathNode[0] = mainNode;
258     for (Int_t i=0; i<=p2volLevel; i++) {
259       pathNode[i+1] = pathNode[i]->GetDaughter(p2nodeInd[i]);
260     };
261     Double_t globalCoord1[3] = {coord1[0], coord1[1], coord1[2]}; 
262     Double_t globalVect1[3]  = {vect1[0], vect1[1], vect1[2]};
263
264     for (Int_t i = commonMotherLevel+1; i <= p2volLevel; i++) {
265       pathNode[i+1]->GetMatrix()->MasterToLocal(globalCoord1, coord1);
266       pathNode[i+1]->GetMatrix()->MasterToLocalVect(globalVect1, vect1);
267       CopyFrom(globalCoord1, coord1);
268       CopyFrom(globalVect1, vect1);
269     };
270   } else {
271     if (! GetCheckPoint(p1, 0, 0, coord1) ) return 0;
272     if (! GetCheckVect(p1, 0, 0, vect1) ) return 0;
273   };
274   
275   //=================================================
276   // Get p2 position in the systeme of p2
277   if (! GetCheckPoint(p2, 0, 0, coord2) ) return 0;
278   if (! GetCheckVect(p2, 0, 0, vect2) ) return 0;
279
280   Double_t cx = (coord1[0]+coord2[0])/2;
281   Double_t cy = (coord1[1]+coord2[1])/2;
282   Double_t cz = (coord1[2]+coord2[2])/2;
283   Double_t dx = coord2[0]-coord1[0];
284   Double_t dy = coord2[1]-coord1[1];
285   Double_t dz = coord2[2]-coord1[2];
286
287   //=================================================
288   // Positionning of the segment between the 2 points
289   if (TMath::Abs(dy)<1e-231) dy = 1e-231;
290   if (TMath::Abs(dz)<1e-231) dz = 1e-231;
291   //Double_t angleRot1 = -TMath::ATan(dx/dy);
292   //Double_t planDiagL = -TMath::Sqrt(dy*dy+dx*dx);
293   //if (dy<0) planDiagL = -planDiagL;
294   //Double_t angleRotDiag = TMath::ATan(planDiagL/dz);
295
296   Double_t angleRot1    = -TMath::ATan2(dx,dy);
297   Double_t planDiagL    =  TMath::Sqrt(dy*dy+dx*dx);
298   Double_t angleRotDiag = -TMath::ATan2(planDiagL,dz);
299   //--- (Calculate rotation of segment on the Z axis)
300   //-- Here I'm trying to calculate the rotation to be applied in
301   //-- order to match as closer as possible this segment and the 
302   //-- previous one. 
303   //-- It seems that some times it doesn't work ...
304   TGeoRotation rotTemp("",angleRot1*TMath::RadToDeg(),
305                        angleRotDiag*TMath::RadToDeg(), rotation);
306   Double_t localX[3] = {0,1,0};
307   Double_t globalX[3];
308   rotTemp.LocalToMasterVect(localX, globalX);
309   CopyFrom(localX, globalX);
310   GetCheckVect(localX, p2Vol, 0, fgkCableMaxNodeLevel+1, globalX);
311   Double_t orthVect[3];
312   GetCheckVect(vect1, p2Vol, 0, fgkCableMaxNodeLevel+1, orthVect);
313 //   Double_t angleRotZ = 0;
314 //   if (p2>1) {
315 //     Double_t orthVectNorm2 = ScalProd(orthVect,orthVect);
316 //     Double_t alpha1 = ScalProd(fPreviousX,orthVect)/orthVectNorm2;
317 //     Double_t alpha2 = ScalProd(globalX,orthVect)/orthVectNorm2;
318 //     Double_t globalX1p[3], globalX2p[3];
319 //     globalX1p[0] = fPreviousX[0] - alpha1*orthVect[0];
320 //     globalX1p[1] = fPreviousX[1] - alpha1*orthVect[1];
321 //     globalX1p[2] = fPreviousX[2] - alpha1*orthVect[2];
322 //     globalX2p[0] = globalX[0] - alpha2*orthVect[0];
323 //     globalX2p[1] = globalX[1] - alpha2*orthVect[1];
324 //     globalX2p[2] = globalX[2] - alpha2*orthVect[2];
325 //     //-- now I'm searching the 3th vect which makes an orthogonal base
326 //     //-- with orthVect and globalX1p ...
327 //     Double_t nulVect[3] = {0,0,0};
328 //     Double_t axis3[3];
329 //     TMath::Normal2Plane(nulVect, orthVect, globalX1p, axis3);
330 //     Double_t globalX1pNorm2 = ScalProd(globalX1p, globalX1p);
331 //     Double_t beta = ScalProd(globalX2p, globalX1p)/globalX1pNorm2;
332 //     Double_t gamma = ScalProd(globalX2p, axis3);
333 //     angleRotZ = (TMath::ATan2(1,0) - TMath::ATan2(beta, gamma))
334 //                 *TMath::RadToDeg();
335 //   };
336   //   cout << "!!!!!!!!!!!!!!!!!!!  angle = " <<angleRotZ << endl;
337   CopyFrom(fPreviousX, globalX);
338   //---
339   Double_t localVect1[3], localVect2[3];
340   TGeoRotation rot("",angleRot1*TMath::RadToDeg(),
341                    angleRotDiag*TMath::RadToDeg(),
342                    rotation);
343 //                 rotation-angleRotZ);
344 // since angleRotZ doesn't always work, I won't use it ...
345
346   rot.MasterToLocalVect(vect1, localVect1);
347   rot.MasterToLocalVect(vect2, localVect2);
348
349   //=================================================
350   // Create the segment and add it to the mother volume
351   TGeoVolume *vCableSegB = CreateSegment(coord1, coord2,
352                                          localVect1, localVect2);
353   TGeoRotation rotArbSeg("", 0, 90, 0);
354   rotArbSeg.MultiplyBy(&rot, kFALSE);
355   TGeoTranslation trans("",cx, cy, cz);
356   TGeoCombiTrans  *combiB = new TGeoCombiTrans(trans, rotArbSeg);
357   p2Vol->AddNode(vCableSegB, p2, combiB);
358   //=================================================;
359
360   if (fDebug) {
361     printf("---\n  Cable segment points : ");
362     printf("%f, %f, %f\n",coord1[0], coord1[1], coord1[2]);
363     printf("%f, %f, %f\n",coord2[0], coord2[1], coord2[2]);
364   };
365
366 //   #include <TGeoSphere.h>
367 //   TGeoMedium *airSDD = gGeoManager->GetMedium("ITS_AIR$");
368 //   TGeoSphere *sphere = new TGeoSphere(0, 0.05);
369 //   TGeoVolume *vSphere = new TGeoVolume("", sphere, airSDD);
370 //   TGeoTranslation *trC = new TGeoTranslation("", cx, cy, cz);
371 //   TGeoTranslation *tr1 = new TGeoTranslation("",coord1[0],
372 //                                           coord1[1],coord1[2]);
373 //   TGeoTranslation *tr2 = new TGeoTranslation("",coord2[0],
374 //                                           coord2[1],coord2[2]);
375 //   p2Vol->AddNode(vSphere, p2*3-2, trC);
376 //   p2Vol->AddNode(vSphere, p2*3-1, tr1);
377 //   p2Vol->AddNode(vSphere, p2*3  , tr2);
378   if (ct) *ct = combiB;
379   return vCableSegB;
380 }
381
382 //________________________________________________________________________
383 TGeoVolume* AliITSv11GeomCableFlat::CreateAndInsertBoxCableSegment(Int_t p2,
384                                                                    Double_t rotation,
385                                                                    TGeoCombiTrans** ct)
386 {
387   // This function is to be use only when the segment has the shape
388   // of a simple box, i.e. the normal vector to its end is perpendicular
389   // to the segment own axis
390 //    Creates a cable segment between points p1 and p2.
391 //    Rotation is the eventual rotation of the flat cable
392 //    along its length axis
393 //
394 // The segment volume is created inside the volume containing point2
395 // Therefore this segment should be defined in this volume only.
396 // I mean here that, if the previous point is in another volume,
397 // it should be just at the border between the 2 volumes. Also the
398 // orientation vector of the previous point should be orthogonal to
399 // the surface between the 2 volumes.
400
401   TGeoNode *mainNode;
402   if (fInitialNode==0) {
403     TObjArray *nodes = gGeoManager->GetListOfNodes();
404     if (nodes->GetEntriesFast()==0) return 0;
405     mainNode = (TGeoNode *) nodes->UncheckedAt(0);
406   } else {
407     mainNode = fInitialNode;
408   };
409
410   Int_t p1 = p2 - 1;
411   TGeoVolume *p2Vol = GetVolume(p2);
412   TGeoVolume *p1Vol = GetVolume(p1);
413
414   ResetCheckDaughter();
415   fCurrentVol = p1Vol;
416   if (! CheckDaughter(mainNode)) {
417     printf("Error::volume containing point is not visible in node tree!\n");
418     return 0;
419   };
420
421   Double_t coord1[3], coord2[3], vect1[3], vect2[3];
422   //=================================================
423   // Get p1 position in the systeme of p2
424   if (p1Vol!=p2Vol) {
425
426     Int_t p1nodeInd[fgkCableMaxNodeLevel]; 
427     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p1nodeInd[i]=fNodeInd[i];
428     Int_t p1volLevel = 0;
429     while (p1nodeInd[p1volLevel]!=-1) p1volLevel++;
430     p1volLevel--;
431
432     ResetCheckDaughter();
433     fCurrentVol = p2Vol;
434     if (! CheckDaughter(mainNode)) {
435       printf("Error::volume containing point is not visible in node tree!\n");
436       return 0;
437     };
438     Int_t p2nodeInd[fgkCableMaxNodeLevel];
439     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p2nodeInd[i]=fNodeInd[i];
440     Int_t commonMotherLevel = 0;
441     while (p1nodeInd[commonMotherLevel]==fNodeInd[commonMotherLevel])
442       commonMotherLevel++;
443     commonMotherLevel--;
444     Int_t p2volLevel = 0;
445     while (fNodeInd[p2volLevel]!=-1) p2volLevel++;
446     p2volLevel--;
447
448     // Get coord and vect of p1 in the common mother reference system
449     if (! GetCheckPoint(p1, 0, p1volLevel-commonMotherLevel, coord1) )
450       return 0;
451     if (! GetCheckVect( p1, 0, p1volLevel-commonMotherLevel, vect1) )
452       return 0;
453
454     // Translate them in the reference system of the volume containing p2    
455     TGeoNode *pathNode[fgkCableMaxNodeLevel];
456     pathNode[0] = mainNode;
457     for (Int_t i=0; i<=p2volLevel; i++) {
458       pathNode[i+1] = pathNode[i]->GetDaughter(p2nodeInd[i]);
459     };
460     Double_t globalCoord1[3] = {coord1[0], coord1[1], coord1[2]}; 
461     Double_t globalVect1[3]  = {vect1[0], vect1[1], vect1[2]};
462
463     for (Int_t i = commonMotherLevel+1; i <= p2volLevel; i++) {
464       pathNode[i+1]->GetMatrix()->MasterToLocal(globalCoord1, coord1);
465       pathNode[i+1]->GetMatrix()->MasterToLocalVect(globalVect1, vect1);
466       CopyFrom(globalCoord1, coord1);
467       CopyFrom(globalVect1, vect1);
468     };
469   } else {
470     if (! GetCheckPoint(p1, 0, 0, coord1) ) return 0;
471     if (! GetCheckVect(p1, 0, 0, vect1) ) return 0;
472   };
473   
474   //=================================================
475   // Get p2 position in the systeme of p2
476   if (! GetCheckPoint(p2, 0, 0, coord2) ) return 0;
477   if (! GetCheckVect(p2, 0, 0, vect2) ) return 0;
478
479   Double_t cx = (coord1[0]+coord2[0])/2;
480   Double_t cy = (coord1[1]+coord2[1])/2;
481   Double_t cz = (coord1[2]+coord2[2])/2;
482   Double_t dx = coord2[0]-coord1[0];
483   Double_t dy = coord2[1]-coord1[1];
484   Double_t dz = coord2[2]-coord1[2];
485
486   //=================================================
487   // Positionning of the segment between the 2 points
488   if (TMath::Abs(dy)<1e-231) dy = 1e-231;
489   if (TMath::Abs(dz)<1e-231) dz = 1e-231;
490   //Double_t angleRot1 = -TMath::ATan(dx/dy);
491   //Double_t planDiagL = -TMath::Sqrt(dy*dy+dx*dx);
492   //if (dy<0) planDiagL = -planDiagL;
493   //Double_t angleRotDiag = TMath::ATan(planDiagL/dz);
494
495   Double_t angleRot1    = -TMath::ATan2(dx,dy);
496   Double_t planDiagL    =  TMath::Sqrt(dy*dy+dx*dx);
497   Double_t angleRotDiag = -TMath::ATan2(planDiagL,dz);
498   //--- (Calculate rotation of segment on the Z axis)
499   //-- Here I'm trying to calculate the rotation to be applied in
500   //-- order to match as closer as possible this segment and the 
501   //-- previous one. 
502   //-- It seems that some times it doesn't work ...
503   TGeoRotation rotTemp("",angleRot1*TMath::RadToDeg(),
504                        angleRotDiag*TMath::RadToDeg(), rotation);
505   Double_t localX[3] = {0,1,0};
506   Double_t globalX[3];
507   rotTemp.LocalToMasterVect(localX, globalX);
508   CopyFrom(localX, globalX);
509   GetCheckVect(localX, p2Vol, 0, fgkCableMaxNodeLevel+1, globalX);
510   Double_t orthVect[3];
511   GetCheckVect(vect1, p2Vol, 0, fgkCableMaxNodeLevel+1, orthVect);
512 //   Double_t angleRotZ = 0;
513 //   if (p2>1) {
514 //     Double_t orthVectNorm2 = ScalProd(orthVect,orthVect);
515 //     Double_t alpha1 = ScalProd(fPreviousX,orthVect)/orthVectNorm2;
516 //     Double_t alpha2 = ScalProd(globalX,orthVect)/orthVectNorm2;
517 //     Double_t globalX1p[3], globalX2p[3];
518 //     globalX1p[0] = fPreviousX[0] - alpha1*orthVect[0];
519 //     globalX1p[1] = fPreviousX[1] - alpha1*orthVect[1];
520 //     globalX1p[2] = fPreviousX[2] - alpha1*orthVect[2];
521 //     globalX2p[0] = globalX[0] - alpha2*orthVect[0];
522 //     globalX2p[1] = globalX[1] - alpha2*orthVect[1];
523 //     globalX2p[2] = globalX[2] - alpha2*orthVect[2];
524 //     //-- now I'm searching the 3th vect which makes an orthogonal base
525 //     //-- with orthVect and globalX1p ...
526 //     Double_t nulVect[3] = {0,0,0};
527 //     Double_t axis3[3];
528 //     TMath::Normal2Plane(nulVect, orthVect, globalX1p, axis3);
529 //     Double_t globalX1pNorm2 = ScalProd(globalX1p, globalX1p);
530 //     Double_t beta = ScalProd(globalX2p, globalX1p)/globalX1pNorm2;
531 //     Double_t gamma = ScalProd(globalX2p, axis3);
532 //     angleRotZ = (TMath::ATan2(1,0) - TMath::ATan2(beta, gamma))
533 //                 *TMath::RadToDeg();
534 //   };
535   CopyFrom(fPreviousX, globalX);
536   //---
537   Double_t localVect1[3], localVect2[3];
538   TGeoRotation rot("",angleRot1*TMath::RadToDeg(),
539                    angleRotDiag*TMath::RadToDeg(),
540                    rotation);
541 //                 rotation-angleRotZ);
542 // since angleRotZ doesn't always work, I won't use it ...
543
544   rot.MasterToLocalVect(vect1, localVect1);
545   rot.MasterToLocalVect(vect2, localVect2);
546
547   //=================================================
548   // Create the segment and add it to the mother volume
549   TGeoVolume *vCableSegB = CreateBoxSegment(coord1, coord2);
550
551   TGeoRotation rotArbSeg("", 0, 90, 0);
552   rotArbSeg.MultiplyBy(&rot, kFALSE);
553   TGeoTranslation trans("",cx, cy, cz);
554   TGeoCombiTrans  *combiB = new TGeoCombiTrans(trans, rotArbSeg);
555   p2Vol->AddNode(vCableSegB, p2, combiB);
556   //=================================================;
557
558   if (fDebug) {
559     printf("---\n  Cable segment points : ");
560     printf("%f, %f, %f\n",coord1[0], coord1[1], coord1[2]);
561     printf("%f, %f, %f\n",coord2[0], coord2[1], coord2[2]);
562   };
563
564   if (ct) *ct = combiB;
565   return vCableSegB;
566 }
567
568 //________________________________________________________________________
569 TGeoVolume* AliITSv11GeomCableFlat::CreateAndInsertCableCylSegment(Int_t p2,
570                                                                    Double_t rotation,
571                                                                    TGeoCombiTrans** ct)
572 {
573   // Create a flat cable segment with a curvature between points p1 and p2.
574   // The radius and position of the curve is defined by the
575   // perpendicular vector of point p2 (the orientation of this vector
576   // and the position of the 2 check points are enough to completely
577   // define the curve)
578   //    Rotation is the eventual rotation of the flat cable
579   //    along its length axis
580   //
581
582   TGeoNode *mainNode;
583   if (fInitialNode==0) {
584     TObjArray *nodes = gGeoManager->GetListOfNodes();
585     if (nodes->GetEntriesFast()==0) return 0;
586     mainNode = (TGeoNode *) nodes->UncheckedAt(0);
587   } else {
588     mainNode = fInitialNode;
589   };
590
591   Int_t p1 = p2 - 1;
592   TGeoVolume *p1Vol = GetVolume(p1);
593   TGeoVolume *p2Vol = GetVolume(p2);
594
595   ResetCheckDaughter();
596   fCurrentVol = p1Vol;
597   if (! CheckDaughter(mainNode)) {
598     printf("Error::volume containing point is not visible in node tree!\n");
599     return 0;
600   };
601
602   Double_t coord1[3], coord2[3], vect1[3], vect2[3];
603   //=================================================
604   // Get p1 position in the systeme of p2
605   if (p1Vol!=p2Vol) {
606
607     Int_t p1nodeInd[fgkCableMaxNodeLevel]; 
608     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p1nodeInd[i]=fNodeInd[i];
609     Int_t p1volLevel = 0;
610     while (p1nodeInd[p1volLevel]!=-1) p1volLevel++;
611     p1volLevel--;
612
613     ResetCheckDaughter();
614     fCurrentVol = p2Vol;
615     if (! CheckDaughter(mainNode)) {
616       printf("Error::volume containing point is not visible in node tree!\n");
617       return 0;
618     };
619     Int_t p2nodeInd[fgkCableMaxNodeLevel];
620     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p2nodeInd[i]=fNodeInd[i];
621     Int_t commonMotherLevel = 0;
622     while (p1nodeInd[commonMotherLevel]==fNodeInd[commonMotherLevel])
623       commonMotherLevel++;
624     commonMotherLevel--;
625     Int_t p2volLevel = 0;
626     while (fNodeInd[p2volLevel]!=-1) p2volLevel++;
627     p2volLevel--;
628
629     // Get coord and vect of p1 in the common mother reference system
630     GetCheckPoint(p1, 0, p1volLevel-commonMotherLevel, coord1);
631     GetCheckVect( p1, 0, p1volLevel-commonMotherLevel, vect1);
632     // Translate them in the reference system of the volume containing p2    
633     TGeoNode *pathNode[fgkCableMaxNodeLevel];
634     pathNode[0] = mainNode;
635     for (Int_t i=0; i<=p2volLevel; i++) {
636       pathNode[i+1] = pathNode[i]->GetDaughter(p2nodeInd[i]);
637     };
638     Double_t globalCoord1[3] = {coord1[0], coord1[1], coord1[2]}; 
639     Double_t globalVect1[3]  = {vect1[0], vect1[1], vect1[2]};
640
641     for (Int_t i = commonMotherLevel+1; i<=p2volLevel; i++) {
642       pathNode[i+1]->GetMatrix()->MasterToLocal(globalCoord1, coord1);
643       pathNode[i+1]->GetMatrix()->MasterToLocalVect(globalVect1, vect1);
644       CopyFrom(globalCoord1, coord1);
645       CopyFrom(globalVect1, vect1);
646     };
647   } else {
648     GetCheckPoint(p1, 0, 0, coord1);
649     GetCheckVect(p1, 0, 0, vect1);
650   };
651   
652   //=================================================
653   // Get p2 position in the systeme of p2
654   GetCheckPoint(p2, 0, 0, coord2);
655   GetCheckVect(p2, 0, 0, vect2);
656
657   Double_t cx = (coord1[0]+coord2[0])/2;
658   Double_t cy = (coord1[1]+coord2[1])/2;
659   Double_t cz = (coord1[2]+coord2[2])/2;
660   Double_t dx = coord2[0]-coord1[0];
661   Double_t dy = coord2[1]-coord1[1];
662   Double_t dz = coord2[2]-coord1[2];
663   Double_t length = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
664
665   //=================================================
666   // Positionning of the segment between the 2 points
667   if ((dy<1e-31)&&(dy>0)) dy = 1e-31;
668   if ((dz<1e-31)&&(dz>0)) dz = 1e-31;
669   if ((dy>-1e-31)&&(dy<0)) dy = -1e-31;
670   if ((dz>-1e-31)&&(dz<0)) dz = -1e-31;
671
672   Double_t angleRot1 = -TMath::ATan2(dx,dy);
673   Double_t planDiagL = TMath::Sqrt(dy*dy+dx*dx);
674   Double_t angleRotDiag = -TMath::ATan2(planDiagL,dz);
675
676   TGeoRotation rotTorusTemp("",angleRot1*TMath::RadToDeg(),
677                             angleRotDiag*TMath::RadToDeg(),0);
678   TGeoRotation rotTorusToZ("",0,90,0);
679   rotTorusTemp.MultiplyBy(&rotTorusToZ, kTRUE);
680   Double_t localVect2[3];
681   rotTorusTemp.MasterToLocalVect(vect2, localVect2);
682   if (localVect2[1]<0) {
683     localVect2[0] = -localVect2[0];
684     localVect2[1] = -localVect2[1];
685     localVect2[2] = -localVect2[2];
686   };
687   Double_t normVect2 = TMath::Sqrt(localVect2[0]*localVect2[0]+
688                                    localVect2[1]*localVect2[1]+
689                                    localVect2[2]*localVect2[2]);
690   Double_t axisX[3] = {1,0,0};
691   Double_t cosangleTorusSeg = (localVect2[0]*axisX[0]+
692                                localVect2[1]*axisX[1]+
693                                localVect2[2]*axisX[2])/normVect2;
694   Double_t angleTorusSeg = TMath::ACos(cosangleTorusSeg)*TMath::RadToDeg();
695   TGeoRotation rotTorus("",angleRot1*TMath::RadToDeg(),
696                         angleRotDiag*TMath::RadToDeg(),
697                         45-angleTorusSeg+rotation);
698                         //180-angleTorusSeg+rotation);
699   rotTorus.MultiplyBy(&rotTorusToZ, kTRUE);
700   rotTorus.MasterToLocalVect(vect2, localVect2);
701   if (localVect2[1]<0) {
702     localVect2[0] = -localVect2[0];
703     localVect2[1] = -localVect2[1];
704     localVect2[2] = -localVect2[2];
705   };
706   normVect2 = TMath::Sqrt(localVect2[0]*localVect2[0]+
707                           localVect2[1]*localVect2[1]+
708                           localVect2[2]*localVect2[2]);
709   Double_t axisY[3] = {0,1,0};
710   Double_t cosPhi = (localVect2[0]*axisY[0]+localVect2[1]*axisY[1]+
711                      localVect2[2]*axisY[2])/normVect2;
712   Double_t torusPhi1 = TMath::ACos(cosPhi);
713   Double_t torusR = (length/2)/TMath::Sin(torusPhi1);
714   torusPhi1 = torusPhi1*TMath::RadToDeg();
715   Double_t perpLength = TMath::Sqrt((torusR-0.5*length)*(torusR+0.5*length));
716   Double_t localTransT[3] = {-perpLength,0,0};
717   Double_t globalTransT[3];
718   rotTorus.LocalToMasterVect(localTransT, globalTransT);
719   TGeoTranslation transTorus("",cx+globalTransT[0],cy+globalTransT[1],
720                              cz+globalTransT[2]);
721
722   TGeoCombiTrans  *combiTorus = new TGeoCombiTrans(transTorus, rotTorus);
723
724   //=================================================
725   // Create the segment and add it to the mother volume
726   TGeoVolume *vCableSegT = CreateCylSegment(torusPhi1, torusR);
727   p2Vol->AddNode(vCableSegT, p2, combiTorus);
728
729   if (fDebug) {
730     printf("---\n  Cable segment points : ");
731     printf("%f, %f, %f\n",coord1[0], coord1[1], coord1[2]);
732     printf("%f, %f, %f\n",coord2[0], coord2[1], coord2[2]);
733   };
734
735   if (ct) *ct = combiTorus;
736   return vCableSegT;
737 }
738
739 //________________________________________________________________________
740 TGeoVolume *AliITSv11GeomCableFlat::CreateSegment( Double_t *coord1,
741                                                    Double_t *coord2,
742                                                    Double_t *localVect1,
743                                                    Double_t *localVect2 )
744 {
745   // Create a segment with arbitrary vertices (general case)
746   //=================================================
747   // Calculate segment "deformation"
748   Double_t dx = coord2[0]-coord1[0];
749   Double_t dy = coord2[1]-coord1[1];
750   Double_t dz = coord2[2]-coord1[2];
751   Double_t length = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
752
753   Double_t cosTheta1 = -1./TMath::Sqrt( 1 + localVect1[0]*localVect1[0]
754                                         /localVect1[2]/localVect1[2] );
755   Double_t cosTheta2 = 1./TMath::Sqrt( 1 + localVect2[0]*localVect2[0]
756                                         /localVect2[2]/localVect2[2] );
757   if (localVect1[2]<0) cosTheta1 = -cosTheta1;
758   if (localVect2[2]<0) cosTheta2 = -cosTheta2;
759
760   Double_t dL1 = 0.5*fWidth*TMath::Tan(TMath::ACos(cosTheta1));
761   Double_t dL2 = 0.5*fWidth*TMath::Tan(TMath::ACos(cosTheta2));
762   if (localVect1[0]<0) dL1 = - dL1;
763   if (localVect2[0]<0) dL2 = - dL2;
764   //---
765   Double_t cosPhi1 = -1./TMath::Sqrt( 1 + localVect1[1]*localVect1[1]
766                                         /localVect1[2]/localVect1[2] );
767   Double_t cosPhi2 = 1./TMath::Sqrt( 1 + localVect2[1]*localVect2[1]
768                                         /localVect2[2]/localVect2[2] );
769   if (localVect1[2]<0) cosPhi1 = -cosPhi1;
770   if (localVect2[2]<0) cosPhi2 = -cosPhi2;
771
772   Double_t tanACosCosPhi1 = TMath::Tan(TMath::ACos(cosPhi1));
773   Double_t tanACosCosPhi2 = TMath::Tan(TMath::ACos(cosPhi2));
774   if (localVect1[1]<0) tanACosCosPhi1 = -tanACosCosPhi1;
775   if (localVect2[1]<0) tanACosCosPhi2 = -tanACosCosPhi2;
776
777   Double_t dl1 = 0.5*fThick*tanACosCosPhi1*0.99999999999999;
778   Double_t dl2 = 0.5*fThick*tanACosCosPhi2*0.99999999999999;
779   // 0.9999999999999 is for correcting problems in TGeo...
780   //=================================================
781   // Create the segment
782   TGeoArb8 *cableSeg = new TGeoArb8(fThick/2);
783   cableSeg->SetVertex( 0, -fWidth/2, -length/2 - dL1 + dl1);
784   cableSeg->SetVertex( 1, -fWidth/2,  length/2 + dL2 - dl2);
785   cableSeg->SetVertex( 2,  fWidth/2,  length/2 - dL2 - dl2);
786   cableSeg->SetVertex( 3,  fWidth/2, -length/2 + dL1 + dl1);
787   cableSeg->SetVertex( 4, -fWidth/2, -length/2 - dL1 - dl1);
788   cableSeg->SetVertex( 5, -fWidth/2,  length/2 + dL2 + dl2);
789   cableSeg->SetVertex( 6,  fWidth/2,  length/2 - dL2 + dl2);
790   cableSeg->SetVertex( 7,  fWidth/2, -length/2 + dL1 - dl1);
791
792   TGeoVolume *vCableSeg = new TGeoVolume(GetName(), cableSeg, fLayMedia[fNlayer-1]);
793   vCableSeg->SetLineColor(fLayColor[fNlayer-1]);
794
795   // add all cable layers but the last
796   for (Int_t iLay=0; iLay<fNlayer-1; iLay++) {
797
798     Double_t dl1Lay = 0.5*fLayThickness[iLay]*tanACosCosPhi1;
799     Double_t dl2Lay = 0.5*fLayThickness[iLay]*tanACosCosPhi2;
800  
801     Double_t ztr = -fThick/2;
802     for (Int_t i=0;i<iLay; i++) ztr+= fLayThickness[i];
803     ztr+= fLayThickness[iLay]/2;
804
805     Double_t dl1LayS = ztr*tanACosCosPhi1;
806     Double_t dl2LayS = ztr*tanACosCosPhi2;
807
808     TGeoArb8 *lay = new TGeoArb8(fLayThickness[iLay]/2);
809     lay->SetVertex( 0, -fWidth/2, -length/2 - dL1 + dl1Lay - dl1LayS);
810     lay->SetVertex( 1, -fWidth/2,  length/2 + dL2 - dl2Lay + dl2LayS);
811     lay->SetVertex( 2,  fWidth/2,  length/2 - dL2 - dl2Lay + dl2LayS);
812     lay->SetVertex( 3,  fWidth/2, -length/2 + dL1 + dl1Lay - dl1LayS);
813     lay->SetVertex( 4, -fWidth/2, -length/2 - dL1 - dl1Lay - dl1LayS);
814     lay->SetVertex( 5, -fWidth/2,  length/2 + dL2 + dl2Lay + dl2LayS);
815     lay->SetVertex( 6,  fWidth/2,  length/2 - dL2 + dl2Lay + dl2LayS);
816     lay->SetVertex( 7,  fWidth/2, -length/2 + dL1 - dl1Lay - dl1LayS);
817     TGeoVolume *vLay = new TGeoVolume("vCableSegLay", lay, fLayMedia[iLay]);
818     vLay->SetLineColor(fLayColor[iLay]);
819     
820     if (fTranslation[iLay]==0)
821       fTranslation[iLay] = new TGeoTranslation(0, 0, ztr);
822     vCableSeg->AddNode(vLay, iLay+1, fTranslation[iLay]);
823   };
824
825   //vCableSeg->SetVisibility(kFALSE);
826   return vCableSeg;
827 }
828
829 //________________________________________________________________________
830 TGeoVolume *AliITSv11GeomCableFlat::CreateCylSegment(Double_t &phi,
831                                                      Double_t &r)
832 {
833   // Create a segment in shape of a cylinder, allows to represent
834   // a folded flat cable
835
836   Double_t phi1 = 360-phi;
837   Double_t phi2 = 360+phi;
838
839   Double_t rMin = r-fThick/2;
840   Double_t rMax = r+fThick/2;
841   //=================================================
842   // Create the segment
843
844   TGeoTubeSeg *cableSeg = new TGeoTubeSeg(rMin, rMax, fWidth/2,
845                                           phi1, phi2);
846   TGeoVolume *vCableSeg = new TGeoVolume(GetName(), cableSeg, fLayMedia[fNlayer-1]);
847   vCableSeg->SetLineColor(fLayColor[fNlayer-1]);
848
849   // add all cable layers but the last
850   for (Int_t iLay=0; iLay<fNlayer-1; iLay++) {
851  
852     Double_t ztr = -fThick/2;
853     for (Int_t i=0;i<iLay; i++) ztr+= fLayThickness[i];
854
855     rMin = r  + ztr;
856     rMax = r  + ztr + fLayThickness[iLay];
857     TGeoTubeSeg *lay = new TGeoTubeSeg(rMin, rMax, fWidth/2,
858                                        phi1, phi2);
859
860     TGeoVolume *vLay = new TGeoVolume("vCableSegLay", lay, fLayMedia[iLay]);
861     vLay->SetLineColor(fLayColor[iLay]);
862     
863     vCableSeg->AddNode(vLay, iLay+1, 0);
864   };
865
866   //vCableSeg->SetVisibility(kFALSE);
867   return vCableSeg;
868 }
869
870 //________________________________________________________________________
871 TGeoVolume *AliITSv11GeomCableFlat::CreateBoxSegment( Double_t *coord1,
872                                                       Double_t *coord2)
873 {
874   // Create a segment for the case it is a simple box
875   //=================================================
876   Double_t dx = coord2[0]-coord1[0];
877   Double_t dy = coord2[1]-coord1[1];
878   Double_t dz = coord2[2]-coord1[2];
879   Double_t length = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
880
881   TGeoBBox *cableSeg = new  TGeoBBox(fWidth/2, length/2, fThick/2);
882   TGeoVolume *vCableSeg = new TGeoVolume(GetName(), cableSeg, fLayMedia[fNlayer-1]);
883   vCableSeg->SetLineColor(fLayColor[fNlayer-1]);
884   // This volume is the cable container. It codes also the material for the
885   // last layer
886
887   // add all cable layers but the last one
888   for (Int_t iLay=0; iLay<fNlayer-1; iLay++) {
889  
890     Double_t ztr = -fThick/2;
891     for (Int_t i=0;i<iLay; i++) ztr+= fLayThickness[i];
892     ztr+= fLayThickness[iLay]/2;
893
894     TGeoBBox *lay = new  TGeoBBox(fWidth/2, length/2, fLayThickness[iLay]/2);
895
896     TGeoVolume *vLay = new TGeoVolume("vCableSegLay", lay, fLayMedia[iLay]);
897     vLay->SetLineColor(fLayColor[iLay]);
898     
899     if (fTranslation[iLay]==0)
900       fTranslation[iLay] = new TGeoTranslation(0, 0, ztr);
901     vCableSeg->AddNode(vLay, iLay+1, fTranslation[iLay]);
902   };
903
904   //vCableSeg->SetVisibility(kFALSE);
905   return vCableSeg;
906 }
907
908 //________________________________________________________________________
909 void AliITSv11GeomCableFlat::SetNLayers(Int_t nLayers) {
910   // Set the number of layers
911   if((nLayers>0) &&(nLayers<=fgkCableMaxLayer)) {
912
913     fNlayer = nLayers;
914     for (Int_t i=0; i<fgkCableMaxLayer ; i++) {
915       fLayThickness[i] = 0;
916       fTranslation[i]  = 0;
917       fLayColor[i]     = 0;
918       fLayMedia[i]     = 0;  
919     }; 
920   };
921 }
922
923 //________________________________________________________________________
924 Int_t AliITSv11GeomCableFlat::SetLayer(Int_t nLayer, Double_t thick,
925                                           TGeoMedium *medium, Int_t color) {
926   // Set the layer number nLayer
927   if ((nLayer<0)||(nLayer>=fNlayer)) {
928     printf("Set wrong layer number of the cable\n");
929     return kFALSE;
930   };
931   if (nLayer>0)
932     if (fLayThickness[nLayer-1]<=0) {
933       printf("AliITSv11GeomCableFlat::SetLayer():"
934              " You must define cable layer %i first !",nLayer-1);
935       return kFALSE;
936     };
937
938   Double_t thickTot = 0;
939   for (Int_t i=0; i<nLayer; i++) thickTot += fLayThickness[i];
940   thickTot += thick;
941   if (thickTot-1e-10>fThick) {
942     printf("Can't add this layer, cable thickness would be higher than total\n");
943     return kFALSE;
944   };
945
946   fLayThickness[nLayer] = thick;
947   fLayMedia[nLayer] = medium;
948   fLayColor[nLayer] = color;
949   fTranslation[nLayer]  = 0;
950   return kTRUE;
951 }