ITS/AliITSv11GeomCableRound.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17
  18 // General Root includes
  19 //#include <Riostream.h>
  20 #include <TMath.h>
  21 #include <TVectorD.h>
  22
  23 // Root Geometry includes
  24 #include <TGeoManager.h>
  25 #include <TGeoVolume.h>
  26 #include <TGeoTube.h>
  27 #include <TGeoMatrix.h>
  28
  29 #include "AliITSv11GeomCableRound.h"
  30
  31 //*************************************************************************
  32 //   Class for round cables
  33 //
  34 // Ludovic Gaudichet                                   gaudichet@to.infn.it
  35 //*************************************************************************
  36
  37 ClassImp(AliITSv11GeomCableRound)
  38
  39 //________________________________________________________________________
  40 AliITSv11GeomCableRound::
  41 AliITSv11GeomCableRound(const char* name, Double_t radius) :
  42 AliITSv11GeomCable(name) {
  43   // Constructor
  44   fRadius = radius;
  45   fNlayer = 0;
  46   for (Int_t i=0; i<fgkCableMaxLayer ; i++) {
  47     fLayThickness[i] = 0;
  48     fLayColor[i] = 0;
  49     fLayMedia[i] = 0;
  50   };
  51   fPhiMin = 0;
  52   fPhiMax = 360;
  53 };
  54
  55 //________________________________________________________________________
  56 AliITSv11GeomCableRound::AliITSv11GeomCableRound(const AliITSv11GeomCableRound &s) :
  57   AliITSv11GeomCable(s),fRadius(s.fRadius),fNlayer(s.fNlayer),fPhiMin(s.fPhiMin),
  58   fPhiMax(s.fPhiMax)
  59 {
  60   //     Copy Constructor
  61   for (Int_t i=0; i<s.fNlayer; i++) {
  62     fLayThickness[i] = s.fLayThickness[i];
  63     fLayMedia[i] = s.fLayMedia[i];
  64     fLayColor[i] = s.fLayColor[i];
  65   }
  66 }
  67
  68 //________________________________________________________________________
  69 AliITSv11GeomCableRound& AliITSv11GeomCableRound::
  70 operator=(const AliITSv11GeomCableRound &s) {
  71   //     Assignment operator
  72   // Not fully inplemented yet !!!
  73
  74   if(&s == this) return *this;
  75   *this = s;
  76   fRadius = s.fRadius;
  77   fPhiMin = s.fPhiMin;
  78   fPhiMax = s.fPhiMax;
  79   fNlayer = s.fNlayer;
  80   for (Int_t i=0; i<s.fNlayer; i++) {
  81     fLayThickness[i] = s.fLayThickness[i];
  82     fLayMedia[i] = s.fLayMedia[i];
  83     fLayColor[i] = s.fLayColor[i];
  84   };
  85   return *this;
  86 }
  87
  88 //________________________________________________________________________
  89 Int_t AliITSv11GeomCableRound::GetPoint( Int_t iCheckPt, Double_t *coord)
  90   const {
  91   // Get check point #iCheckPt
  92   TVectorD *coordVector =(TVectorD *)fPointArray.UncheckedAt(2*iCheckPt);
  93   CopyFrom(coord, coordVector->GetMatrixArray());
  94   return kTRUE;
  95 };
  96
  97 //________________________________________________________________________
  98 Int_t AliITSv11GeomCableRound::GetVect( Int_t iCheckPt, Double_t *coord)
  99   const {
 100   //
 101   // Get vector transverse to the section at point #iCheckPt
 102   //
 103
 104   TVectorD *coordVector =(TVectorD *)fPointArray.UncheckedAt(2*iCheckPt+1);
 105   CopyFrom(coord, coordVector->GetMatrixArray());
 106   return kTRUE;
 107 };
 108 //________________________________________________________________________
 109 void AliITSv11GeomCableRound::AddCheckPoint( TGeoVolume *vol, Int_t iCheckPt,
 110                                                Double_t *coord, Double_t *orthVect)
 111 {
 112   //
 113   // Add point #iCheckPt and its transverse vector. Point is added at (i) in
 114   // fPointArray and the vector is added at (i+1)
 115   //
 116
 117
 118   if (iCheckPt>=fVolumeArray.GetEntriesFast()) {
 119     fVolumeArray.AddLast(vol);
 120     TVectorD *point = new TVectorD(3,coord);
 121     TVectorD *vect  = new TVectorD(3,orthVect);
 122     fPointArray.AddLast(point);
 123     fPointArray.AddLast(vect);
 124
 125   } else if ((iCheckPt >= 0)&&(iCheckPt < fVolumeArray.GetEntriesFast())) {
 126     fVolumeArray.AddAt(vol, iCheckPt);
 127     TVectorD *point = new TVectorD(3,coord);
 128     TVectorD *vect  = new TVectorD(3,orthVect);
 129     fPointArray.AddAt(point, iCheckPt*2  );
 130     fPointArray.AddAt(vect,  iCheckPt*2+1);
 131   };
 132 };
 133
 134 //________________________________________________________________________
 135 void AliITSv11GeomCableRound::PrintCheckPoints() const {
 136   // Print all check points
 137
 138   printf("  ---\n  Printing all check points of the round cable\n");
 139   for (Int_t i = 0; i<fVolumeArray.GetEntriesFast(); i++) {
 140     TVectorD *coordVector = (TVectorD *)fPointArray.UncheckedAt(i*2);
 141     //TVectorD *vectVector = (TVectorD *)fPointArray.UncheckedAt(i*2+1);
 142     Double_t coord[3];
 143     CopyFrom(coord, coordVector->GetMatrixArray());
 144
 145     printf("   ( %.2f, %.2f, %.2f )\n", coord[0], coord[1], coord[2]);
 146   };
 147
 148 };
 149
 150 //________________________________________________________________________
 151 Int_t AliITSv11GeomCableRound::CreateAndInsertCableSegment(Int_t p2)
 152 {
 153 //    Creates a cable segment between points p1 and p2.
 154 //    Rotation is the eventual rotation of the flat cable
 155 //    along its length axis
 156 //
 157 // The segment volume is created inside the volume containing point2
 158 // Therefore this segment should be defined in this volume only.
 159 // I mean here that, if the previous point is in another volume,
 160 // it should be just at the border between the 2 volumes. Also the
 161 // orientation vector of the previous point should be othogonal to
 162 // the surface between the 2 volumes.
 163
 164   TGeoNode *mainNode;
 165   if (fInitialNode==0) {
 166     TObjArray *nodes = gGeoManager->GetListOfNodes();
 167     if (nodes->GetEntriesFast()==0) return kFALSE;
 168     mainNode = (TGeoNode *) nodes->UncheckedAt(0);
 169   } else {
 170     mainNode = fInitialNode;
 171   };
 172
 173   Int_t p1 = p2 - 1;
 174   TGeoVolume *p1Vol = GetVolume(p1);
 175   TGeoVolume *p2Vol = GetVolume(p2);
 176
 177   ResetCheckDaughter();
 178   fCurrentVol = p1Vol;
 179   if (! CheckDaughter(mainNode)) {
 180     printf("Error::volume containing point is not visible in node tree!\n");
 181     return kFALSE;
 182   };
 183
 184   Double_t coord1[3], coord2[3], vect1[3], vect2[3];
 185   //=================================================
 186   // Get p1 position in the systeme of p2
 187   if (p1Vol!=p2Vol) {
 188
 189     Int_t p1nodeInd[fgkCableMaxNodeLevel];
 190     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p1nodeInd[i]=fNodeInd[i];
 191     Int_t p1volLevel = 0;
 192     while (p1nodeInd[p1volLevel]!=-1) p1volLevel++;
 193     p1volLevel--;
 194
 195     ResetCheckDaughter();
 196     fCurrentVol = p2Vol;
 197     if (! CheckDaughter(mainNode)) {
 198       printf("Error::volume containing point is not visible in node tree!\n");
 199       return kFALSE;
 200     };
 201     Int_t p2nodeInd[fgkCableMaxNodeLevel];
 202     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p2nodeInd[i]=fNodeInd[i];
 203     Int_t commonMotherLevel = 0;
 204     while (p1nodeInd[commonMotherLevel]==fNodeInd[commonMotherLevel])
 205       commonMotherLevel++;
 206     commonMotherLevel--;
 207     Int_t p2volLevel = 0;
 208     while (fNodeInd[p2volLevel]!=-1) p2volLevel++;
 209     p2volLevel--;
 210
 211     // Get coord and vect of p1 in the common mother reference system
 212     GetCheckPoint(p1, 0, p1volLevel-commonMotherLevel, coord1);
 213     GetCheckVect( p1, 0, p1volLevel-commonMotherLevel, vect1);
 214     // Translate them in the reference system of the volume containing p2
 215     TGeoNode *pathNode[fgkCableMaxNodeLevel];
 216     pathNode[0] = mainNode;
 217     for (Int_t i=0; i<=p2volLevel; i++) {
 218       pathNode[i+1] = pathNode[i]->GetDaughter(p2nodeInd[i]);
 219     };
 220     Double_t globalCoord1[3] = {coord1[0], coord1[1], coord1[2]};
 221     Double_t globalVect1[3]  = {vect1[0], vect1[1], vect1[2]};
 222
 223     for (Int_t i = commonMotherLevel+1; i<=p2volLevel; i++) {
 224       pathNode[i+1]->GetMatrix()->MasterToLocal(globalCoord1, coord1);
 225       pathNode[i+1]->GetMatrix()->MasterToLocalVect(globalVect1, vect1);
 226       CopyFrom(globalCoord1, coord1);
 227       CopyFrom(globalVect1, vect1);
 228     };
 229   } else {
 230     GetCheckPoint(p1, 0, 0, coord1);
 231     GetCheckVect(p1, 0, 0, vect1);
 232   };
 233
 234   //=================================================
 235   // Get p2 position in the systeme of p2
 236   GetCheckPoint(p2, 0, 0, coord2);
 237   GetCheckVect(p2, 0, 0, vect2);
 238
 239   Double_t cx = (coord1[0]+coord2[0])/2;
 240   Double_t cy = (coord1[1]+coord2[1])/2;
 241   Double_t cz = (coord1[2]+coord2[2])/2;
 242   Double_t dx = coord2[0]-coord1[0];
 243   Double_t dy = coord2[1]-coord1[1];
 244   Double_t dz = coord2[2]-coord1[2];
 245
 246   //=================================================
 247   // Positionning of the segment between the 2 points
 248   if ((dy<1e-31)&&(dy>0)) dy = 1e-31;
 249   if ((dz<1e-31)&&(dz>0)) dz = 1e-31;
 250   if ((dy>-1e-31)&&(dy<0)) dy = -1e-31;
 251   if ((dz>-1e-31)&&(dz<0)) dz = -1e-31;
 252
 253   Double_t angleRot1 = -TMath::ATan2(dx,dy);
 254   Double_t planDiagL = TMath::Sqrt(dy*dy+dx*dx);
 255   Double_t angleRotDiag = -TMath::ATan2(planDiagL,dz);
 256   TGeoRotation *rot = new TGeoRotation("",angleRot1*TMath::RadToDeg(),
 257                                        angleRotDiag*TMath::RadToDeg(),
 258                                        0);
 259   Double_t localVect1[3], localVect2[3];
 260   rot->MasterToLocalVect(vect1, localVect1);
 261   rot->MasterToLocalVect(vect2, localVect2);
 262   TGeoTranslation *trans = new TGeoTranslation("",cx, cy, cz);
 263
 264   //=================================================
 265   // Create the segment and add it to the mother volume
 266   TGeoVolume *vCableSeg = CreateSegment(coord1, coord2,
 267                                         localVect1, localVect2);
 268
 269   TGeoCombiTrans  *combi = new TGeoCombiTrans(*trans, *rot);
 270   p2Vol->AddNode(vCableSeg, p2, combi);
 271   //=================================================
 272   delete rot;
 273   delete trans;
 274
 275   if (fDebug) {
 276     printf("---\n  Cable segment points : ");
 277     printf("%f, %f, %f\n",coord1[0], coord1[1], coord1[2]);
 278     printf("%f, %f, %f\n",coord2[0], coord2[1], coord2[2]);
 279   };
 280 //   #include <TGeoSphere.h>
 281 //   TGeoMedium *airSDD = gGeoManager->GetMedium("ITSsddAir");
 282 //   TGeoSphere *sphere = new TGeoSphere(0, 0.15);
 283 //   TGeoVolume *vSphere = new TGeoVolume("", sphere, airSDD);
 284 //   TGeoTranslation *trC = new TGeoTranslation("", cx, cy, cz);
 285 //   TGeoTranslation *tr1 = new TGeoTranslation("",coord1[0],
 286 //                                           coord1[1],coord1[2]);
 287 //   TGeoTranslation *tr2 = new TGeoTranslation("",coord2[0],
 288 //                                           coord2[1],coord2[2]);
 289 //   p2Vol->AddNode(vSphere, p2*3-2, trC);
 290 //   p2Vol->AddNode(vSphere, p2*3-1, tr1);
 291 //   p2Vol->AddNode(vSphere, p2*3  , tr2);
 292
 293   return kTRUE;
 294 };
 295
 296 //________________________________________________________________________
 297 TGeoVolume *AliITSv11GeomCableRound::CreateSegment( Double_t *coord1,
 298                                                       Double_t *coord2,
 299                                                       Double_t *localVect1,
 300                                                       Double_t *localVect2 )
 301 {
 302
 303   //=================================================
 304   // Calculate segment "deformation"
 305   Double_t dx = coord2[0]-coord1[0];
 306   Double_t dy = coord2[1]-coord1[1];
 307   Double_t dz = coord2[2]-coord1[2];
 308   Double_t length = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
 309
 310   // normal vectors have to point outside the TGeoCtub :
 311   if (-localVect1[2]<0) {
 312     localVect1[0] = -localVect1[0];
 313     localVect1[1] = -localVect1[1];
 314     localVect1[2] = -localVect1[2];
 315   };
 316   if (localVect2[2]<0) {
 317     localVect2[0] = -localVect2[0];
 318     localVect2[1] = -localVect2[1];
 319     localVect2[2] = -localVect2[2];
 320   };
 321   //=================================================
 322   // Create the segment
 323   TGeoCtub *cableSeg = new TGeoCtub(0., fRadius, length/2, fPhiMin, fPhiMax,
 324                                     localVect1[0],localVect1[1],localVect1[2],
 325                                     localVect2[0],localVect2[1],localVect2[2]);
 326
 327   TGeoMedium *airSDD = gGeoManager->GetMedium("ITSair");
 328   TGeoVolume *vCableSeg = new TGeoVolume(GetName(), cableSeg, airSDD);
 329
 330   // add all cable layers
 331   Double_t layThickness[100+1];                        // 100 layers max !!!
 332   layThickness[0] = 0;
 333   for (Int_t iLay=0; iLay<fNlayer; iLay++) {
 334
 335     layThickness[iLay+1] = fLayThickness[iLay]+layThickness[iLay];
 336     TGeoCtub *lay = new TGeoCtub(layThickness[iLay], layThickness[iLay+1],
 337                                  length/2, fPhiMin, fPhiMax,
 338                                  localVect1[0],localVect1[1],localVect1[2],
 339                                  localVect2[0],localVect2[1],localVect2[2]);
 340
 341     TGeoVolume *vLay = new TGeoVolume("vCableSegLay", lay, fLayMedia[iLay]);
 342     vLay->SetLineColor(fLayColor[iLay]);
 343     vCableSeg->AddNode(vLay, iLay+1, 0);
 344   };
 345
 346   vCableSeg->SetVisibility(kFALSE);
 347   return vCableSeg;
 348 };
 349
 350
 351 //________________________________________________________________________
 352 void AliITSv11GeomCableRound::SetNLayers(Int_t nLayers) {
 353   // Set the total number of layers
 354   if((nLayers>0) &&(nLayers<=fgkCableMaxLayer)) {
 355     fNlayer = nLayers;
 356     for (Int_t i = 0; i<fNlayer; i++) {
 357       fLayThickness[i] = 0;
 358       fLayMedia[i] = 0;
 359     };
 360   };
 361 };
 362
 363 //________________________________________________________________________
 364 Int_t AliITSv11GeomCableRound::SetLayer(Int_t nLayer, Double_t thick,
 365                                            TGeoMedium *medium, Int_t color) {
 366   // Set layer #nLayer
 367   if ((nLayer<0)||(nLayer>=fNlayer)) {
 368     printf("Set wrong layer number of the cable\n");
 369     return kFALSE;
 370   };
 371   if (nLayer>0)
 372     if (fLayThickness[nLayer-1]<=0) {
 373       printf("You must define cable layer %i first !",nLayer-1);
 374       return kFALSE;
 375     };
 376
 377   Double_t thickTot = 0;
 378   for (Int_t i=0; i<nLayer; i++) thickTot += fLayThickness[i];
 379   thickTot += thick;
 380   if (thickTot-1e-10>fRadius) {
 381     printf("Can't add this layer, cable thickness would be higher than total\n");
 382     return kFALSE;
 383   };
 384
 385   fLayThickness[nLayer] = thick;
 386   fLayMedia[nLayer] = medium;
 387   fLayColor[nLayer] = color;
 388
 389   return kTRUE;
 390 };