New geometry: SDD, cables and update on V11 (L. Gaudichet)
[u/mrichter/AliRoot.git] / ITS / AliITSv11GeomCableFlat.cxx
CommitLineData
b7943f00 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
40ClassImp(AliITSv11GeomCableFlat)
41
42//________________________________________________________________________
43AliITSv11GeomCableFlat::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//________________________________________________________________________
58AliITSv11GeomCableFlat::
59AliITSv11GeomCableFlat(const char* name, Double_t width, Double_t thick) :
60AliITSv11GeomCable(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//________________________________________________________________________
74AliITSv11GeomCableFlat::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//________________________________________________________________________
87AliITSv11GeomCableFlat& AliITSv11GeomCableFlat::
88operator=(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//________________________________________________________________________
107Int_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//________________________________________________________________________
121Int_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//________________________________________________________________________
136void 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//________________________________________________________________________
166void 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//________________________________________________________________________
177Int_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//________________________________________________________________________
373TGeoVolume *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//________________________________________________________________________
464void 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//________________________________________________________________________
479Int_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};