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