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