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