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