d7a2e62ce8fa177e5a8ce49e9ac5c6c1bb9b630f
[u/mrichter/AliRoot.git] / FMD / AliFMDGeometry.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 /* $Id$ */
16 /** @file    AliFMDGeometry.cxx
17     @author  Christian Holm Christensen <cholm@nbi.dk>
18     @date    Mon Mar 27 12:40:37 2006
19     @brief   Geometry mananger for the FMD
20 */
21 //____________________________________________________________________
22 //                                                                          
23 // Forward Multiplicity Detector based on Silicon wafers. 
24 //
25 // This class is a singleton that handles the geometry parameters of
26 // the FMD detectors.  
27 //                                                       
28 // The actual code is done by various separate classes.   Below is
29 // diagram showing the relationship between the various FMD classes
30 // that handles the geometry 
31 //
32 //                               +------------+ 
33 //                            +- | AliFMDRing |
34 //                         2  |  +------------+
35 //      +----------------+<>--+        |                                
36 //      | AliFMDGeometry |             ^                        
37 //      +----------------+<>--+        V 1..2                           
38 //                         3  | +----------------+              
39 //                            +-| AliFMDDetector |              
40 //                              +----------------+              
41 //                                     ^
42 //                                     |
43 //                       +-------------+-------------+
44 //                       |             |             |        
45 //                  +---------+   +---------+   +---------+
46 //                  | AliFMD1 |   | AliFMD2 |   | AliFMD3 |
47 //                  +---------+   +---------+   +---------+
48 //      
49 //
50 // *  AliFMDRing 
51 //    This class contains all stuff needed to do with a ring.  It's
52 //    used by the AliFMDDetector objects to instantise inner and
53 //    outer rings.  The AliFMDRing objects are shared by the
54 //    AliFMDDetector objects, and owned by the AliFMDv1 object. 
55 //
56 // *  AliFMD1, AliFMD2, and AliFMD3 
57 //    These are specialisation of AliFMDDetector, that contains the
58 //    particularities of each of the sub-detector system.  It is
59 //    envisioned that the classes should also define the support
60 //    volumes and material for each of the detectors.                          
61 //                                                                          
62 //
63 #include "AliFMDGeometry.h"     // ALIFMDGEOMETRY_H
64 #include "AliFMDRing.h"         // ALIFMDRING_H
65 #include "AliFMD1.h"            // ALIFMD1_H
66 #include "AliFMD2.h"            // ALIFMD2_H
67 #include "AliFMD3.h"            // ALIFMD2_H
68 #include "AliRecPoint.h"        // ALIRECPOINT_H
69 #include "AliFMDDebug.h"                   // ALILOG_H
70 #include <TVector3.h>           // ROOT_TVector3
71 // #include <TMatrix.h>            // ROOT_TMatrix
72 // #include <TParticle.h>          // ROOT_TParticle
73 #include <Riostream.h>
74 #include "AliFMDGeometryBuilder.h"
75 // #include <TArrayI.h>
76 #include <TGeoManager.h>
77 #include <TGeoVolume.h>
78 #include <TGeoNode.h>
79 #include <TMath.h>
80 static Int_t FindNodeDepth(const char* name, const char* volname);
81
82
83 //====================================================================
84 ClassImp(AliFMDGeometry)
85 #if 0
86   ; // This is here to keep Emacs for indenting the next line
87 #endif
88
89 //____________________________________________________________________
90 AliFMDGeometry* AliFMDGeometry::fgInstance = 0;
91
92 //____________________________________________________________________
93 AliFMDGeometry* 
94 AliFMDGeometry::Instance() 
95 {
96   // Return (newly created) singleton instance 
97   if (!fgInstance) fgInstance = new AliFMDGeometry;
98   return fgInstance;
99 }
100
101 //____________________________________________________________________
102 AliFMDGeometry::AliFMDGeometry() 
103   : AliGeometry("FMD", "Forward multiplicity"), 
104     fIsInitialized(kFALSE), 
105     fInner(0),
106     fOuter(0),
107     fFMD1(0),
108     fFMD2(0),
109     fFMD3(0),
110     fUseFMD1(kFALSE),
111     fUseFMD2(kFALSE),
112     fUseFMD3(kFALSE),
113     fIsInitTrans(kFALSE),
114     fBuilder(0),
115     fDetectorOff(0),
116     fModuleOff(0),  
117     fRingOff(0),
118     fSectorOff(0),
119     fActive(2),
120     fDetailed(kFALSE),       
121     fUseAssembly(kFALSE)
122 {
123   // PROTECTED
124   // Default constructor 
125   fUseFMD1     = kTRUE;
126   fUseFMD2     = kTRUE;
127   fUseFMD3     = kTRUE;  
128   fDetailed    = kTRUE;
129   fUseAssembly = kTRUE;
130   fInner = new AliFMDRing('I');
131   fOuter = new AliFMDRing('O');
132   fFMD1  = new AliFMD1(fInner);
133   fFMD2  = new AliFMD2(fInner, fOuter);
134   fFMD3  = new AliFMD3(fInner, fOuter);
135   fIsInitialized = kFALSE;
136   fActive.Reset(-1);
137 }
138
139 //____________________________________________________________________
140 AliFMDGeometry::AliFMDGeometry(const AliFMDGeometry& other) 
141   : AliGeometry(other),
142     fIsInitialized(other.fIsInitialized),
143     fInner(other.fInner), 
144     fOuter(other.fOuter), 
145     fFMD1(other.fFMD1), 
146     fFMD2(other.fFMD2), 
147     fFMD3(other.fFMD3), 
148     fUseFMD1(other.fUseFMD1), 
149     fUseFMD2(other.fUseFMD2), 
150     fUseFMD3(other.fUseFMD3), 
151     fIsInitTrans(other.fIsInitTrans),
152     fBuilder(other.fBuilder),
153     fDetectorOff(other.fDetectorOff),
154     fModuleOff(other.fModuleOff),  
155     fRingOff(other.fRingOff),
156     fSectorOff(other.fSectorOff),
157     fActive(other.fActive),
158     fDetailed(other.fDetailed),
159     fUseAssembly(other.fUseAssembly)
160 {
161   // PROTECTED
162   // Copy constructor
163 }
164
165
166
167 //____________________________________________________________________
168 AliFMDGeometry&
169 AliFMDGeometry::operator=(const AliFMDGeometry& other) 
170 {
171   // PROTECTED
172   // Assignment operator 
173   fUseFMD1              = other.fUseFMD1; 
174   fUseFMD2              = other.fUseFMD2; 
175   fUseFMD3              = other.fUseFMD3; 
176   fFMD1                 = other.fFMD1; 
177   fFMD2                 = other.fFMD2; 
178   fFMD3                 = other.fFMD3; 
179   fInner                = other.fInner; 
180   fOuter                = other.fOuter; 
181   fIsInitialized        = other.fIsInitialized;
182   return *this;
183 }
184
185 //____________________________________________________________________
186 void
187 AliFMDGeometry::Init()
188 {
189   // Initialize the the singleton if not done so already 
190   if (fIsInitialized) return;
191   fInner->Init();
192   fOuter->Init();
193   fFMD1->Init();
194   fFMD2->Init();
195   fFMD3->Init();
196 }
197
198 //____________________________________________________________________
199 void
200 AliFMDGeometry::InitTransformations(Bool_t force)
201 {
202   // Find all local <-> global transforms 
203   if (force) fIsInitTrans = kFALSE;
204   if (fIsInitTrans) return; 
205   if (!gGeoManager) {
206     AliError("No TGeoManager defined");
207     return;
208   }
209   AliFMDDebug(1, ("Initialising transforms for FMD geometry"));
210   if (fFMD1) fFMD1->InitTransformations();
211   if (fFMD2) fFMD2->InitTransformations();
212   if (fFMD3) fFMD3->InitTransformations();
213   fIsInitTrans = kTRUE;
214 }
215
216 //____________________________________________________________________
217 void
218 AliFMDGeometry::Build()
219 {
220   // Build the geometry 
221   if (!fBuilder) fBuilder = new AliFMDGeometryBuilder(fDetailed);
222   fBuilder->SetDetailed(fDetailed);
223   fBuilder->UseAssembly(fUseAssembly);
224   fBuilder->Exec();
225 }
226
227 //____________________________________________________________________
228 void
229 AliFMDGeometry::SetActive(Int_t* active, Int_t n) 
230 {
231   // Set active volumes 
232   fActive.Set(n);
233   for (Int_t i = 0; i < n; i++) { 
234     AliFMDDebug(1, ("Active vol id # %d: %d", i, active[i]));
235     fActive[i] = active[i];
236   }
237 }
238
239 //____________________________________________________________________
240 void
241 AliFMDGeometry::AddActive(Int_t active)
242 {
243   // Add an active volume 
244   Int_t n = fActive.fN;
245   fActive.Set(n+1);
246   fActive[n] = active;
247 }
248
249 //____________________________________________________________________
250 Bool_t
251 AliFMDGeometry::IsActive(Int_t vol) const
252 {
253   // Check if a volume is active 
254   for (Int_t i = 0; i < fActive.fN; i++) 
255     if (fActive[i] == vol) return kTRUE;
256   return kFALSE;
257 }
258   
259 //____________________________________________________________________
260 AliFMDDetector*
261 AliFMDGeometry::GetDetector(Int_t i) const
262 {
263   // Get the ith detector.   i should be one of 1, 2, or 3.  If an
264   // invalid value is passed, 0 (NULL) is returned. 
265   switch (i) {
266   case 1: return fUseFMD1 ? static_cast<AliFMDDetector*>(fFMD1) : 0;
267   case 2: return fUseFMD2 ? static_cast<AliFMDDetector*>(fFMD2) : 0;
268   case 3: return fUseFMD3 ? static_cast<AliFMDDetector*>(fFMD3) : 0;
269   }
270   return 0;
271 }
272
273 //____________________________________________________________________
274 AliFMDRing*
275 AliFMDGeometry::GetRing(Char_t i) const
276 {
277   // Get the ith ring.  i should be one of 'I' or 'O' (case
278   // insensitive).  If an invalid parameter is passed, 0 (NULL) is
279   // returned. 
280   switch (i) {
281   case 'I':
282   case 'i': return fInner;
283   case 'O':
284   case 'o': return fOuter;
285   }
286   return 0;
287 }
288
289 //____________________________________________________________________
290 void
291 AliFMDGeometry::Enable(Int_t i)
292 {
293   // Enable the ith detector.  i should be one of 1, 2, or 3
294   switch (i) {
295   case 1: fUseFMD1 = kTRUE; break;
296   case 2: fUseFMD2 = kTRUE; break;
297   case 3: fUseFMD3 = kTRUE; break;
298   }
299 }
300
301 //____________________________________________________________________
302 void
303 AliFMDGeometry::Disable(Int_t i)
304 {
305   // Disable the ith detector.  i should be one of 1, 2, or 3
306   switch (i) {
307   case 1: fUseFMD1 = kFALSE; break;
308   case 2: fUseFMD2 = kFALSE; break;
309   case 3: fUseFMD3 = kFALSE; break;
310   }
311 }
312
313 //____________________________________________________________________
314 void
315 AliFMDGeometry::Detector2XYZ(UShort_t  detector, 
316                              Char_t    ring, 
317                              UShort_t  sector, 
318                              UShort_t  strip, 
319                              Double_t& x, 
320                              Double_t& y, 
321                              Double_t& z) const
322 {
323   // Translate detector coordinates (detector, ring, sector, strip) to
324   // spatial coordinates (x, y, z) in the master reference frame of
325   // ALICE. 
326   AliFMDDetector* det = GetDetector(detector);
327   if (!det) { 
328     AliWarning(Form("Unknown detector %d", detector));
329     return;
330   }
331   det->Detector2XYZ(ring, sector, strip, x, y, z);
332 }
333
334 //____________________________________________________________________
335 Bool_t
336 AliFMDGeometry::XYZ2Detector(Double_t  x, 
337                              Double_t  y, 
338                              Double_t  z,
339                              UShort_t& detector, 
340                              Char_t&   ring, 
341                              UShort_t& sector, 
342                              UShort_t& strip) const
343 {
344   // Translate spatial coordinates (x,y,z) in the master reference frame of
345   // ALICE to the detector coordinates (detector, ring, sector,
346   // strip).  Note, that if this method is to be used in
347   // reconstruction or the like, then the input z-coordinate should be
348   // corrected for the events interactions points z-coordinate, like 
349   // geom->XYZ2Detector(x,y,z-ipz,d,r,s,t);
350   AliFMDDetector* det = 0;
351   detector = 0;
352   for (int i = 1; i <= 3; i++) {
353     det = GetDetector(i);
354     if (!det) continue;
355     if (det->XYZ2Detector(x, y, z, ring, sector, strip)) {
356       detector = det->GetId();
357       return kTRUE;
358     }
359   }
360   return kFALSE;
361 }
362
363 //____________________________________________________________________
364 Bool_t
365 AliFMDGeometry::XYZ2REtaPhiTheta(Double_t  x,   Double_t y, 
366                                  Double_t  z, 
367                                  Double_t& r,   Double_t& eta, 
368                                  Double_t& phi, Double_t& theta)
369 {
370   
371   // Service function to convert Cartisean XYZ to r, eta, phi, and theta.   
372   // 
373   // Note, that the z input should be corrected for the vertex location 
374   // if needed.
375   //
376   //     x      Cartisean X coordinate
377   //     y      Cartisean Y coordinate 
378   //     z      Cartisean Z coordinate 
379   //     r      On return, the radius
380   //     eta    On return, the pseudo-rapidity
381   //     phi    On return, the azimuthal angle
382   //     theta  On return, the polar angle;
383   // 
384   // Return:
385   //     kFALSE in case of problems. 
386
387   if (x == 0 && y == 0 && z == 0) return kFALSE;
388   
389   // Correct for vertex offset. 
390   phi   =  TMath::ATan2(y, x);
391   r     =  TMath::Sqrt(y * y + x * x);
392   theta =  TMath::ATan2(r, z);
393   eta   = -TMath::Log(TMath::Tan(theta / 2));
394
395   return kTRUE;
396 }
397
398
399 //____________________________________________________________________
400 void
401 AliFMDGeometry::GetGlobal(const AliRecPoint* p, 
402                           TVector3& pos, 
403                           TMatrixF& /* mat */) const 
404 {
405   // Get the global coordinates cooresponding to the reconstructed
406   // point p.  The coordiates is returned in the 3-vector pos passed
407   // to this member function.  The matrix mat is used for rotations. 
408   GetGlobal(p, pos);
409 }
410
411 //____________________________________________________________________
412 void
413 AliFMDGeometry::GetGlobal(const AliRecPoint* p, TVector3& pos) const 
414 {
415   // Get the global coordinates cooresponding to the reconstructed
416   // point p.  The coordiates is returned in the 3-vector pos passed
417   // to this member function. Note, as AliRecPoint only has places for
418   // 3 indicies, it is assumed that the ring hit is an inner ring -
419   // which obviously needn't be the case. This makes the member
420   // function pretty darn useless. 
421   // FIXME: Implement this function to work with outer rings too. 
422   Double_t x, y, z;
423   TVector3 local;
424   p->GetLocalPosition(local);
425   UShort_t detector = UShort_t(local.X());
426   UShort_t sector   = UShort_t(local.Y());
427   UShort_t strip    = UShort_t(local.Z());
428   Detector2XYZ(detector, 'I', sector, strip, x, y, z);
429   pos.SetXYZ(x, y, z);
430 }
431
432 //____________________________________________________________________
433 Bool_t
434 AliFMDGeometry::Impact(const TParticle* /* particle */) const 
435
436   // Return true, if the particle will hit the active detector
437   // elements, and false if not.  Should be used for fast
438   // simulations.  Note, that the function currently return false
439   // always.  
440   // FIXME: Implement this function. 
441   return kFALSE; 
442 }
443
444 //____________________________________________________________________  
445 void  
446 AliFMDGeometry::SetAlignableVolumes() const
447 {
448   // Declare alignable volumes
449   for (Int_t d = 1; d <= 3; d++) 
450     if (GetDetector(d)) GetDetector(d)->SetAlignableVolumes();
451 }
452
453
454 //____________________________________________________________________  
455 void  
456 AliFMDGeometry::ExtractGeomInfo()
457 {
458   // Check the volume depth of some nodes, get the active volume
459   // numbers, and so forth. 
460   // 
461   // TODO: Here, we should actually also get the parameters of the
462   // shapes, like the verticies of the polygon shape that makes up the
463   // silicon sensor, the strip pitch, the ring radii, the z-positions,
464   // and so on - that is, all the geometric information we need for
465   // futher processing, such as simulation, digitization,
466   // reconstruction, etc. 
467   Int_t detectorDepth = FindNodeDepth("F1MT_1", "ALIC");
468   Int_t ringDepth     = FindNodeDepth(Form("FITV_%d", int('I')), "ALIC");
469   Int_t moduleDepth   = FindNodeDepth("FIBH_0", "ALIC");
470   Int_t sectorDepth   = FindNodeDepth("FISC_1", "ALIC");
471   fActive.Set(0);
472   fActive.Reset(-1);
473   AliFMDDebug(1, ("Geometry depths:\n"
474                    "   Sector:     %d\n"
475                    "   Module:     %d\n"
476                    "   Ring:       %d\n"
477                    "   Detector:   %d", 
478                    sectorDepth, moduleDepth, ringDepth, detectorDepth));
479   if (sectorDepth < 0 && moduleDepth < 0) {
480     fDetailed    = kFALSE;
481     fSectorOff   = -1;
482     fModuleOff   = -1;
483     fRingOff     = 0;
484     fDetectorOff = (ringDepth - detectorDepth);
485     TGeoVolume* actiVol = gGeoManager->GetVolume("FIAC");
486     TGeoVolume* actoVol = gGeoManager->GetVolume("FOAC");
487     if (actiVol) AddActive(actiVol->GetNumber());
488     if (actiVol) AddActive(actoVol->GetNumber());
489   }
490   else if (sectorDepth < 0) {
491     fDetailed    = kFALSE;
492     fSectorOff   = -1;
493     fModuleOff   = 1;
494     fRingOff     = (moduleDepth - ringDepth) + 1;
495     fDetectorOff = (moduleDepth - detectorDepth) + 1;
496     TGeoVolume* modiVol = gGeoManager->GetVolume("FIMO");
497     TGeoVolume* modoVol = gGeoManager->GetVolume("FOMO");
498     if (modiVol) AddActive(modiVol->GetNumber());
499     if (modoVol) AddActive(modoVol->GetNumber());
500   }
501   else {
502     Int_t stripDepth    = FindNodeDepth("FIST_1", "ALIC");
503     fDetailed    = kTRUE;
504     fSectorOff   = (stripDepth - sectorDepth);
505     fModuleOff   = (moduleDepth >= 0 ? (stripDepth - moduleDepth) : -1);
506     fRingOff     = (stripDepth - ringDepth);
507     fDetectorOff = (stripDepth - detectorDepth );
508     TGeoVolume* striVol = gGeoManager->GetVolume("FIST");
509     TGeoVolume* stroVol = gGeoManager->GetVolume("FOST");
510     if (striVol) AddActive(striVol->GetNumber());
511     if (stroVol) AddActive(stroVol->GetNumber());
512   }    
513   AliFMDDebug(1, ("Geometry offsets:\n"
514                    "   Sector:     %d\n"
515                    "   Module:     %d\n"
516                    "   Ring:       %d\n"
517                    "   Detector:   %d", 
518                    fSectorOff, fModuleOff, fRingOff, fDetectorOff));
519 }
520
521   
522 //____________________________________________________________________  
523 static Int_t 
524 CheckNodes(TGeoNode* node, const char* name, Int_t& lvl)
525 {
526   // If there's no node here. 
527   if (!node) return -1;
528   // Check if it this one 
529   TString sname(name);
530   if (sname == node->GetName()) return lvl;
531
532   // Check if the node is an immediate daugther 
533   TObjArray* nodes = node->GetNodes();
534   if (!nodes) return -1;
535   // Increase the level, and search immediate sub nodes. 
536   lvl++;
537   TGeoNode*  found = static_cast<TGeoNode*>(nodes->FindObject(name));
538   if (found) return lvl;
539
540   // Check the sub node, if any of their sub-nodes match.
541   for (Int_t i = 0; i < nodes->GetEntries(); i++) {
542     TGeoNode* sub = static_cast<TGeoNode*>(nodes->At(i));
543     if (!sub) continue;
544     // Recurive check 
545     if (CheckNodes(sub, name, lvl) >= 0) return lvl;
546   }
547   // If not found, decrease the level 
548   lvl--;
549   return -1;
550 }
551 //____________________________________________________________________  
552 Int_t 
553 FindNodeDepth(const char* name, const char* volname) 
554 {
555   // Find the depth of a node 
556   TGeoVolume* vol  = gGeoManager->GetVolume(volname);
557   if (!vol) {
558     std::cerr << "No top volume defined" << std::endl;
559     return -1;
560   }
561
562   TGeoIterator next(vol);
563   TGeoNode*    node = 0;
564   TString      sName(name);
565   while ((node = next())) { 
566     if (sName == node->GetName()) { 
567       //std::cout << "Found node " << node->GetName() << " at level " 
568       //                << next.GetLevel() << std::endl;
569       return next.GetLevel();
570     }
571   }
572   return -1;
573 #if 0
574   TObjArray* nodes = vol->GetNodes();
575   if (!nodes) { 
576     std::cerr << "No nodes in top volume" << std::endl;
577     return -1;
578   }
579   TIter next(nodes);
580   TGeoNode* node = 0;
581   Int_t lvl = 0;
582   while ((node = static_cast<TGeoNode*>(next()))) 
583     if (CheckNodes(node, name, lvl) >= 0) return lvl;
584   return -1;
585 #endif
586 }
587
588 //____________________________________________________________________
589 //
590 // EOF
591 //