New geometry: SDD, cables and update on V11 (L. Gaudichet)
[u/mrichter/AliRoot.git] / ITS / AliITSv11GeomCableFlat.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
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 <TGeoMatrix.h>
35 #include <TGeoNode.h>
36
37 #include "AliITSv11GeomCableFlat.h"
38
39
40 ClassImp(AliITSv11GeomCableFlat)
41
42 //________________________________________________________________________
43 AliITSv11GeomCableFlat::AliITSv11GeomCableFlat() : AliITSv11GeomCable()
44 {
45   // constructor
46   fWidth = 0;
47   fThick = 0;
48   fNlayer = 0;
49   for (Int_t i=0; i<fgkCableMaxLayer ; i++) {
50     fLayThickness[i] = 0;
51     fTranslation[i]  = 0;
52     fLayColor[i]     = 0;
53     fLayMedia[i]     = 0;  
54  };
55 };
56
57 //________________________________________________________________________
58 AliITSv11GeomCableFlat::
59 AliITSv11GeomCableFlat(const char* name, Double_t width, Double_t thick) :
60 AliITSv11GeomCable(name) {
61   // standard constructor
62   fWidth = width;
63   fThick = thick;
64   fNlayer = 0;
65   for (Int_t i=0; i<fgkCableMaxLayer ; i++) {
66     fLayThickness[i] = 0;
67     fTranslation[i]  = 0;
68     fLayColor[i]     = 0;
69     fLayMedia[i]     = 0;  
70   }; 
71 };
72
73 //________________________________________________________________________
74 AliITSv11GeomCableFlat::AliITSv11GeomCableFlat(const AliITSv11GeomCableFlat &s) :
75   AliITSv11GeomCable(s),fWidth(s.fWidth),fThick(s.fThick),fNlayer(s.fNlayer)
76 {
77   //     Copy Constructor 
78   for (Int_t i=0; i<s.fNlayer; i++) {
79     fLayThickness[i] = s.fLayThickness[i];
80     fTranslation[i] = s.fTranslation[i];
81     fLayMedia[i] = s.fLayMedia[i];
82     fLayColor[i] = s.fLayColor[i];
83   }
84 }
85
86 //________________________________________________________________________
87 AliITSv11GeomCableFlat& AliITSv11GeomCableFlat::
88 operator=(const AliITSv11GeomCableFlat &s) {
89   //     Assignment operator
90   // Not fully inplemented yet !!!
91
92   if(&s == this) return *this;
93   *this = s;
94   fWidth = s.fWidth;
95   fThick = s.fThick;
96   fNlayer = s.fNlayer;
97   for (Int_t i=0; i<s.fNlayer; i++) {
98     fLayThickness[i] = s.fLayThickness[i];
99     fTranslation[i] = s.fTranslation[i];
100     fLayMedia[i] = s.fLayMedia[i];
101     fLayColor[i] = s.fLayColor[i];
102   };
103   return *this;
104 }
105
106 //________________________________________________________________________
107 Int_t AliITSv11GeomCableFlat::GetPoint( Int_t iCheckPt, Double_t *coord)
108   const {
109   // Get the correct point #iCheckPt
110   TVectorD *coordVector =(TVectorD *)fPointArray.At(2*iCheckPt);
111   if (coordVector) {
112     CopyFrom(coord, coordVector->GetMatrixArray());
113     return kTRUE;
114   } else {
115     return kFALSE;
116   };
117 };
118
119
120 //________________________________________________________________________
121 Int_t AliITSv11GeomCableFlat::GetVect( Int_t iCheckPt, Double_t *coord)
122   const {
123   // Get the correct vect corresponding to point #iCheckPt
124
125   TVectorD *coordVector =(TVectorD *)fPointArray.At(2*iCheckPt+1);
126   if (coordVector) {
127     CopyFrom(coord, coordVector->GetMatrixArray());
128     return kTRUE;
129   } else {
130     return kFALSE;
131   };
132 };
133
134
135 //________________________________________________________________________
136 void AliITSv11GeomCableFlat::AddCheckPoint( TGeoVolume *vol, Int_t iCheckPt,
137                                             Double_t *coord, Double_t *orthVect)
138 {
139   //
140   // Add a check point. In the fPointArray, the point is at i and its vector
141   // is at i+1.
142   //
143
144 //   if (iCheckPt>=fVolumeArray.GetEntriesFast()) {
145 //     fVolumeArray.AddLast(vol);
146 //     TVectorD *point = new TVectorD(3,coord);
147 //     TVectorD *vect  = new TVectorD(3,orthVect);
148 //     fPointArray.AddLast(point);
149 //     fPointArray.AddLast(vect);
150
151 //   } else if ((iCheckPt >= 0)&&(iCheckPt < fVolumeArray.GetEntriesFast())) {
152 //     fVolumeArray.AddAt(vol, iCheckPt);
153 //     TVectorD *point = new TVectorD(3,coord);
154 //     TVectorD *vect  = new TVectorD(3,orthVect);
155 //     fPointArray.AddAt(point, iCheckPt*2  );
156 //     fPointArray.AddAt(vect,  iCheckPt*2+1);
157 //   };
158   fVolumeArray.AddAtAndExpand(vol, iCheckPt);
159   TVectorD *point = new TVectorD(3,coord);
160   TVectorD *vect  = new TVectorD(3,orthVect);
161   fPointArray.AddAtAndExpand(point, iCheckPt*2  );
162   fPointArray.AddAtAndExpand(vect,  iCheckPt*2+1);
163 };
164
165 //________________________________________________________________________
166 void AliITSv11GeomCableFlat::PrintCheckPoints() const {
167   // print all check points of the cable
168   printf("  ---\n  Printing all check points of the flat cable\n");
169   for (Int_t i = 0; i<fVolumeArray.GetEntriesFast(); i++) {
170      Double_t coord[3];
171      if (GetPoint( i, coord))
172        printf("   ( %.2f, %.2f, %.2f )\n", coord[0], coord[1], coord[2]);
173   };
174 };
175
176 //________________________________________________________________________
177 Int_t AliITSv11GeomCableFlat::CreateAndInsertCableSegment(Int_t p2,
178                                                           Double_t rotation)
179 {
180 //    Creates a cable segment between points p1 and p2.
181 //    Rotation is the eventual rotation of the flat cable
182 //    along its length axis
183 //
184 // The segment volume is created inside the volume containing point2
185 // Therefore this segment should be defined in this volume only.
186 // I mean here that, if the previous point is in another volume,
187 // it should be just at the border between the 2 volumes. Also the
188 // orientation vector of the previous point should be orthogonal to
189 // the surface between the 2 volumes.
190
191   TGeoNode *mainNode;
192   if (fInitialNode==0) {
193     TObjArray *nodes = gGeoManager->GetListOfNodes();
194     if (nodes->GetEntriesFast()==0) return kFALSE;
195     mainNode = (TGeoNode *) nodes->UncheckedAt(0);
196   } else {
197     mainNode = fInitialNode;
198   };
199
200   Int_t p1 = p2 - 1;
201   TGeoVolume *p2Vol = GetVolume(p2);
202   TGeoVolume *p1Vol = GetVolume(p1);
203
204   ResetCheckDaughter();
205   fCurrentVol = p1Vol;
206   if (! CheckDaughter(mainNode)) {
207     printf("Error::volume containing point is not visible in node tree!\n");
208     return kFALSE;
209   };
210
211   Double_t coord1[3], coord2[3], vect1[3], vect2[3];
212   //=================================================
213   // Get p1 position in the systeme of p2
214   if (p1Vol!=p2Vol) {
215
216     Int_t p1nodeInd[fgkCableMaxNodeLevel]; 
217     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p1nodeInd[i]=fNodeInd[i];
218     Int_t p1volLevel = 0;
219     while (p1nodeInd[p1volLevel]!=-1) p1volLevel++;
220     p1volLevel--;
221
222     ResetCheckDaughter();
223     fCurrentVol = p2Vol;
224     if (! CheckDaughter(mainNode)) {
225       printf("Error::volume containing point is not visible in node tree!\n");
226       return kFALSE;
227     };
228     Int_t p2nodeInd[fgkCableMaxNodeLevel];
229     for (Int_t i=0; i<fgkCableMaxNodeLevel; i++) p2nodeInd[i]=fNodeInd[i];
230     Int_t commonMotherLevel = 0;
231     while (p1nodeInd[commonMotherLevel]==fNodeInd[commonMotherLevel])
232       commonMotherLevel++;
233     commonMotherLevel--;
234     Int_t p2volLevel = 0;
235     while (fNodeInd[p2volLevel]!=-1) p2volLevel++;
236     p2volLevel--;
237
238     // Get coord and vect of p1 in the common mother reference system
239     if (! GetCheckPoint(p1, 0, p1volLevel-commonMotherLevel, coord1) )
240       return kFALSE;
241     if (! GetCheckVect( p1, 0, p1volLevel-commonMotherLevel, vect1) )
242       return kFALSE;
243
244     // Translate them in the reference system of the volume containing p2    
245     TGeoNode *pathNode[fgkCableMaxNodeLevel];
246     pathNode[0] = mainNode;
247     for (Int_t i=0; i<=p2volLevel; i++) {
248       pathNode[i+1] = pathNode[i]->GetDaughter(p2nodeInd[i]);
249     };
250     Double_t globalCoord1[3] = {coord1[0], coord1[1], coord1[2]}; 
251     Double_t globalVect1[3]  = {vect1[0], vect1[1], vect1[2]};
252
253     for (Int_t i = commonMotherLevel+1; i <= p2volLevel; i++) {
254       pathNode[i+1]->GetMatrix()->MasterToLocal(globalCoord1, coord1);
255       pathNode[i+1]->GetMatrix()->MasterToLocalVect(globalVect1, vect1);
256       CopyFrom(globalCoord1, coord1);
257       CopyFrom(globalVect1, vect1);
258     };
259   } else {
260     if (! GetCheckPoint(p1, 0, 0, coord1) ) return kFALSE;
261     if (! GetCheckVect(p1, 0, 0, vect1) ) return kFALSE;
262   };
263   
264   //=================================================
265   // Get p2 position in the systeme of p2
266   if (! GetCheckPoint(p2, 0, 0, coord2) ) return kFALSE;
267   if (! GetCheckVect(p2, 0, 0, vect2) ) return kFALSE;
268
269   Double_t cx = (coord1[0]+coord2[0])/2;
270   Double_t cy = (coord1[1]+coord2[1])/2;
271   Double_t cz = (coord1[2]+coord2[2])/2;
272   Double_t dx = coord2[0]-coord1[0];
273   Double_t dy = coord2[1]-coord1[1];
274   Double_t dz = coord2[2]-coord1[2];
275
276   //=================================================
277   // Positionning of the segment between the 2 points
278   if (TMath::Abs(dy)<1e-231) dy = 1e-231;
279   if (TMath::Abs(dz)<1e-231) dz = 1e-231;
280   //Double_t angleRot1 = -TMath::ATan(dx/dy);
281   //Double_t planDiagL = -TMath::Sqrt(dy*dy+dx*dx);
282   //if (dy<0) planDiagL = -planDiagL;
283   //Double_t angleRotDiag = TMath::ATan(planDiagL/dz);
284
285   Double_t angleRot1    = -TMath::ATan2(dx,dy);
286   Double_t planDiagL    =  TMath::Sqrt(dy*dy+dx*dx);
287   Double_t angleRotDiag = -TMath::ATan2(planDiagL,dz);
288   //--- (Calculate rotation of segment on the Z axis)
289   //-- Here I'm trying to calculate the rotation to be applied in
290   //-- order to match as closer as possible this segment and the 
291   //-- previous one. 
292   //-- It seems that some times it doesn't work ...
293   Double_t angleRotZ = 0;
294   TGeoRotation rotTemp("",angleRot1*TMath::RadToDeg(),
295                        angleRotDiag*TMath::RadToDeg(), rotation);
296   Double_t localX[3] = {0,1,0};
297   Double_t globalX[3];
298   rotTemp.LocalToMasterVect(localX, globalX);
299   CopyFrom(localX, globalX);
300   GetCheckVect(localX, p2Vol, 0, fgkCableMaxNodeLevel+1, globalX);
301   Double_t orthVect[3];
302   GetCheckVect(vect1, p2Vol, 0, fgkCableMaxNodeLevel+1, orthVect);
303   if (p2>1) {
304     Double_t orthVectNorm2 = ScalProd(orthVect,orthVect);
305     Double_t alpha1 = ScalProd(fPreviousX,orthVect)/orthVectNorm2;
306     Double_t alpha2 = ScalProd(globalX,orthVect)/orthVectNorm2;
307     Double_t globalX1p[3], globalX2p[3];
308     globalX1p[0] = fPreviousX[0] - alpha1*orthVect[0];
309     globalX1p[1] = fPreviousX[1] - alpha1*orthVect[1];
310     globalX1p[2] = fPreviousX[2] - alpha1*orthVect[2];
311     globalX2p[0] = globalX[0] - alpha2*orthVect[0];
312     globalX2p[1] = globalX[1] - alpha2*orthVect[1];
313     globalX2p[2] = globalX[2] - alpha2*orthVect[2];
314     //-- now I'm searching the 3th vect which makes an orthogonal base
315     //-- with orthVect and globalX1p ...
316     Double_t nulVect[3] = {0,0,0};
317     Double_t axis3[3];
318     TMath::Normal2Plane(nulVect, orthVect, globalX1p, axis3);
319     Double_t globalX1pNorm2 = ScalProd(globalX1p, globalX1p);
320     Double_t beta = ScalProd(globalX2p, globalX1p)/globalX1pNorm2;
321     Double_t gamma = ScalProd(globalX2p, axis3);
322     angleRotZ = (TMath::ATan2(1,0) - TMath::ATan2(beta, gamma))
323                 *TMath::RadToDeg();
324   };
325   //   cout << "!!!!!!!!!!!!!!!!!!!  angle = " <<angleRotZ << endl;
326   CopyFrom(fPreviousX, globalX);
327   //---
328   Double_t localVect1[3], localVect2[3];
329   TGeoRotation rot("",angleRot1*TMath::RadToDeg(),
330                    angleRotDiag*TMath::RadToDeg(),
331                    rotation);
332 //                 rotation-angleRotZ);
333 // since angleRotZ doesn't always work, I won't use it ...
334
335   rot.MasterToLocalVect(vect1, localVect1);
336   rot.MasterToLocalVect(vect2, localVect2);
337
338   //=================================================
339   // Create the segment and add it to the mother volume
340   TGeoVolume *vCableSegB = CreateSegment(coord1, coord2,
341                                          localVect1, localVect2);
342
343   TGeoRotation rotArbSeg("", 0, 90, 0);
344   rotArbSeg.MultiplyBy(&rot, kFALSE);
345   TGeoTranslation trans("",cx, cy, cz);
346   TGeoCombiTrans  *combiB = new TGeoCombiTrans(trans, rotArbSeg);
347   p2Vol->AddNode(vCableSegB, p2, combiB);
348   //=================================================;
349
350   if (fDebug) {
351     printf("---\n  Cable segment points : ");
352     printf("%f, %f, %f\n",coord1[0], coord1[1], coord1[2]);
353     printf("%f, %f, %f\n",coord2[0], coord2[1], coord2[2]);
354   };
355
356 //   #include <TGeoSphere.h>
357 //   TGeoMedium *airSDD = gGeoManager->GetMedium("ITSsddAir");
358 //   TGeoSphere *sphere = new TGeoSphere(0, 0.05);
359 //   TGeoVolume *vSphere = new TGeoVolume("", sphere, airSDD);
360 //   TGeoTranslation *trC = new TGeoTranslation("", cx, cy, cz);
361 //   TGeoTranslation *tr1 = new TGeoTranslation("",coord1[0],
362 //                                           coord1[1],coord1[2]);
363 //   TGeoTranslation *tr2 = new TGeoTranslation("",coord2[0],
364 //                                           coord2[1],coord2[2]);
365 //   p2Vol->AddNode(vSphere, p2*3-2, trC);
366 //   p2Vol->AddNode(vSphere, p2*3-1, tr1);
367 //   p2Vol->AddNode(vSphere, p2*3  , tr2);
368
369   return kTRUE;
370 };
371
372 //________________________________________________________________________
373 TGeoVolume *AliITSv11GeomCableFlat::CreateSegment( Double_t *coord1,
374                                                       Double_t *coord2,
375                                                       Double_t *localVect1,
376                                                       Double_t *localVect2 )
377 {
378
379   //=================================================
380   // Calculate segment "deformation"
381   Double_t dx = coord2[0]-coord1[0];
382   Double_t dy = coord2[1]-coord1[1];
383   Double_t dz = coord2[2]-coord1[2];
384   Double_t length = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
385
386   Double_t cosTheta1 = -1./TMath::Sqrt( 1 + localVect1[0]*localVect1[0]
387                                         /localVect1[2]/localVect1[2] );
388   Double_t cosTheta2 = 1./TMath::Sqrt( 1 + localVect2[0]*localVect2[0]
389                                         /localVect2[2]/localVect2[2] );
390   if (localVect1[2]<0) cosTheta1 = -cosTheta1;
391   if (localVect2[2]<0) cosTheta2 = -cosTheta2;
392
393   Double_t dL1 = 0.5*fWidth*TMath::Tan(TMath::ACos(cosTheta1));
394   Double_t dL2 = 0.5*fWidth*TMath::Tan(TMath::ACos(cosTheta2));
395   if (localVect1[0]<0) dL1 = - dL1;
396   if (localVect2[0]<0) dL2 = - dL2;
397   //---
398   Double_t cosPhi1 = -1./TMath::Sqrt( 1 + localVect1[1]*localVect1[1]
399                                         /localVect1[2]/localVect1[2] );
400   Double_t cosPhi2 = 1./TMath::Sqrt( 1 + localVect2[1]*localVect2[1]
401                                         /localVect2[2]/localVect2[2] );
402   if (localVect1[2]<0) cosPhi1 = -cosPhi1;
403   if (localVect2[2]<0) cosPhi2 = -cosPhi2;
404
405   Double_t tanACosCosPhi1 = TMath::Tan(TMath::ACos(cosPhi1));
406   Double_t tanACosCosPhi2 = TMath::Tan(TMath::ACos(cosPhi2));
407   if (localVect1[1]<0) tanACosCosPhi1 = -tanACosCosPhi1;
408   if (localVect2[1]<0) tanACosCosPhi2 = -tanACosCosPhi2;
409
410   Double_t dl1 = 0.5*fThick*tanACosCosPhi1;
411   Double_t dl2 = 0.5*fThick*tanACosCosPhi2;
412
413   //=================================================
414   // Create the segment
415   TGeoArb8 *cableSeg = new TGeoArb8(fThick/2);
416   cableSeg->SetVertex( 0, -fWidth/2, -length/2 - dL1 + dl1);
417   cableSeg->SetVertex( 1,  fWidth/2, -length/2 + dL1 + dl1);
418   cableSeg->SetVertex( 2,  fWidth/2,  length/2 - dL2 - dl2);
419   cableSeg->SetVertex( 3, -fWidth/2,  length/2 + dL2 - dl2);
420   cableSeg->SetVertex( 4, -fWidth/2, -length/2 - dL1 - dl1);
421   cableSeg->SetVertex( 5,  fWidth/2, -length/2 + dL1 - dl1);
422   cableSeg->SetVertex( 6,  fWidth/2,  length/2 - dL2 + dl2);
423   cableSeg->SetVertex( 7, -fWidth/2,  length/2 + dL2 + dl2);
424
425   TGeoMedium *airSDD = gGeoManager->GetMedium("ITSair");
426   TGeoVolume *vCableSeg = new TGeoVolume(GetName(), cableSeg, airSDD);
427
428   // add all cable layers
429   for (Int_t iLay=0; iLay<fNlayer; iLay++) {
430
431     Double_t dl1Lay = 0.5*fLayThickness[iLay]*tanACosCosPhi1;
432     Double_t dl2Lay = 0.5*fLayThickness[iLay]*tanACosCosPhi2;
433  
434     Double_t ztr = -fThick/2;
435     for (Int_t i=0;i<iLay; i++) ztr+= fLayThickness[i];
436     ztr+= fLayThickness[iLay]/2;
437
438     Double_t dl1LayS = ztr*tanACosCosPhi1;
439     Double_t dl2LayS = ztr*tanACosCosPhi2;
440
441     TGeoArb8 *lay = new TGeoArb8(fLayThickness[iLay]/2);
442     lay->SetVertex( 0, -fWidth/2, -length/2 - dL1 + dl1Lay - dl1LayS);
443     lay->SetVertex( 1,  fWidth/2, -length/2 + dL1 + dl1Lay - dl1LayS);
444     lay->SetVertex( 2,  fWidth/2,  length/2 - dL2 - dl2Lay + dl2LayS);
445     lay->SetVertex( 3, -fWidth/2,  length/2 + dL2 - dl2Lay + dl2LayS);
446     lay->SetVertex( 4, -fWidth/2, -length/2 - dL1 - dl1Lay - dl1LayS);
447     lay->SetVertex( 5,  fWidth/2, -length/2 + dL1 - dl1Lay - dl1LayS);
448     lay->SetVertex( 6,  fWidth/2,  length/2 - dL2 + dl2Lay + dl2LayS);
449     lay->SetVertex( 7, -fWidth/2,  length/2 + dL2 + dl2Lay + dl2LayS);
450     TGeoVolume *vLay = new TGeoVolume("vCableSegLay", lay, fLayMedia[iLay]);
451     vLay->SetLineColor(fLayColor[iLay]);
452     
453     if (fTranslation[iLay]==0)
454       fTranslation[iLay] = new TGeoTranslation(0, 0, ztr);
455     vCableSeg->AddNode(vLay, iLay+1, fTranslation[iLay]);
456   };
457
458   vCableSeg->SetVisibility(kFALSE);
459   return vCableSeg;
460 };
461
462
463 //________________________________________________________________________
464 void AliITSv11GeomCableFlat::SetNLayers(Int_t nLayers) {
465   // Set the number of layers
466   if((nLayers>0) &&(nLayers<=fgkCableMaxLayer)) {
467
468     fNlayer = nLayers;
469     for (Int_t i=0; i<fgkCableMaxLayer ; i++) {
470       fLayThickness[i] = 0;
471       fTranslation[i]  = 0;
472       fLayColor[i]     = 0;
473       fLayMedia[i]     = 0;  
474     }; 
475   };
476 };
477
478 //________________________________________________________________________
479 Int_t AliITSv11GeomCableFlat::SetLayer(Int_t nLayer, Double_t thick,
480                                           TGeoMedium *medium, Int_t color) {
481   // Set the layer number nLayer
482   if ((nLayer<0)||(nLayer>=fNlayer)) {
483     printf("Set wrong layer number of the cable\n");
484     return kFALSE;
485   };
486   if (nLayer>0)
487     if (fLayThickness[nLayer-1]<=0) {
488       printf("AliITSv11GeomCableFlat::SetLayer():"
489              " You must define cable layer %i first !",nLayer-1);
490       return kFALSE;
491     };
492
493   Double_t thickTot = 0;
494   for (Int_t i=0; i<nLayer; i++) thickTot += fLayThickness[i];
495   thickTot += thick;
496   if (thickTot-1e-10>fThick) {
497     printf("Can't add this layer, cable thickness would be higher than total\n");
498     return kFALSE;
499   };
500
501   fLayThickness[nLayer] = thick;
502   fLayMedia[nLayer] = medium;
503   fLayColor[nLayer] = color;
504   fTranslation[nLayer]  = 0;
505   return kTRUE;
506 };