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 }