Fixes from cholm to enable reading of geometry from a file
[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 static Int_t FindNodeDepth(const char* name, const char* volname);
80
81
82 //====================================================================
83 ClassImp(AliFMDGeometry)
84 #if 0
85   ; // This is here to keep Emacs for indenting the next line
86 #endif
87
88 //____________________________________________________________________
89 AliFMDGeometry* AliFMDGeometry::fgInstance = 0;
90
91 //____________________________________________________________________
92 AliFMDGeometry* 
93 AliFMDGeometry::Instance() 
94 {
95   // Return (newly created) singleton instance 
96   if (!fgInstance) fgInstance = new AliFMDGeometry;
97   return fgInstance;
98 }
99
100 //____________________________________________________________________
101 AliFMDGeometry::AliFMDGeometry() 
102   : AliGeometry("FMD", "Forward multiplicity"), 
103     fIsInitialized(kFALSE), 
104     fInner(0),
105     fOuter(0),
106     fFMD1(0),
107     fFMD2(0),
108     fFMD3(0),
109     fUseFMD1(kFALSE),
110     fUseFMD2(kFALSE),
111     fUseFMD3(kFALSE),
112     fIsInitTrans(kFALSE),
113     fBuilder(0),
114     fDetectorOff(0),
115     fModuleOff(0),  
116     fRingOff(0),
117     fSectorOff(0),
118     fActive(2),
119     fDetailed(kFALSE),       
120     fUseAssembly(kFALSE)
121 {
122   // PROTECTED
123   // Default constructor 
124   fUseFMD1     = kTRUE;
125   fUseFMD2     = kTRUE;
126   fUseFMD3     = kTRUE;  
127   fDetailed    = kTRUE;
128   fUseAssembly = kTRUE;
129   fInner = new AliFMDRing('I');
130   fOuter = new AliFMDRing('O');
131   fFMD1  = new AliFMD1(fInner);
132   fFMD2  = new AliFMD2(fInner, fOuter);
133   fFMD3  = new AliFMD3(fInner, fOuter);
134   fIsInitialized = kFALSE;
135   fActive.Reset(-1);
136 }
137
138 //____________________________________________________________________
139 AliFMDGeometry::AliFMDGeometry(const AliFMDGeometry& other) 
140   : AliGeometry(other),
141     fIsInitialized(other.fIsInitialized),
142     fInner(other.fInner), 
143     fOuter(other.fOuter), 
144     fFMD1(other.fFMD1), 
145     fFMD2(other.fFMD2), 
146     fFMD3(other.fFMD3), 
147     fUseFMD1(other.fUseFMD1), 
148     fUseFMD2(other.fUseFMD2), 
149     fUseFMD3(other.fUseFMD3), 
150     fIsInitTrans(other.fIsInitTrans),
151     fBuilder(other.fBuilder),
152     fDetectorOff(other.fDetectorOff),
153     fModuleOff(other.fModuleOff),  
154     fRingOff(other.fRingOff),
155     fSectorOff(other.fSectorOff),
156     fActive(other.fActive),
157     fDetailed(other.fDetailed),
158     fUseAssembly(other.fUseAssembly)
159 {
160   // PROTECTED
161   // Copy constructor
162 }
163
164
165
166 //____________________________________________________________________
167 AliFMDGeometry&
168 AliFMDGeometry::operator=(const AliFMDGeometry& other) 
169 {
170   // PROTECTED
171   // Assignment operator 
172   fUseFMD1              = other.fUseFMD1; 
173   fUseFMD2              = other.fUseFMD2; 
174   fUseFMD3              = other.fUseFMD3; 
175   fFMD1                 = other.fFMD1; 
176   fFMD2                 = other.fFMD2; 
177   fFMD3                 = other.fFMD3; 
178   fInner                = other.fInner; 
179   fOuter                = other.fOuter; 
180   fIsInitialized        = other.fIsInitialized;
181   return *this;
182 }
183
184 //____________________________________________________________________
185 void
186 AliFMDGeometry::Init()
187 {
188   // Initialize the the singleton if not done so already 
189   if (fIsInitialized) return;
190   fInner->Init();
191   fOuter->Init();
192   fFMD1->Init();
193   fFMD2->Init();
194   fFMD3->Init();
195 }
196
197 //____________________________________________________________________
198 void
199 AliFMDGeometry::InitTransformations(Bool_t force)
200 {
201   // Find all local <-> global transforms 
202   if (force) fIsInitTrans = kFALSE;
203   if (fIsInitTrans) return; 
204   if (!gGeoManager) {
205     AliError("No TGeoManager defined");
206     return;
207   }
208   AliFMDDebug(1, ("Initialising transforms for FMD geometry"));
209   if (fFMD1) fFMD1->InitTransformations();
210   if (fFMD2) fFMD2->InitTransformations();
211   if (fFMD3) fFMD3->InitTransformations();
212   fIsInitTrans = kTRUE;
213 }
214
215 //____________________________________________________________________
216 void
217 AliFMDGeometry::Build()
218 {
219   // Build the geometry 
220   if (!fBuilder) fBuilder = new AliFMDGeometryBuilder(fDetailed);
221   fBuilder->SetDetailed(fDetailed);
222   fBuilder->UseAssembly(fUseAssembly);
223   fBuilder->Exec();
224 }
225
226 //____________________________________________________________________
227 void
228 AliFMDGeometry::SetActive(Int_t* active, Int_t n) 
229 {
230   // Set active volumes 
231   fActive.Set(n);
232   for (Int_t i = 0; i < n; i++) { 
233     AliFMDDebug(1, ("Active vol id # %d: %d", i, active[i]));
234     fActive[i] = active[i];
235   }
236 }
237
238 //____________________________________________________________________
239 void
240 AliFMDGeometry::AddActive(Int_t active)
241 {
242   // Add an active volume 
243   Int_t n = fActive.fN;
244   fActive.Set(n+1);
245   fActive[n] = active;
246 }
247
248 //____________________________________________________________________
249 Bool_t
250 AliFMDGeometry::IsActive(Int_t vol) const
251 {
252   // Check if a volume is active 
253   for (Int_t i = 0; i < fActive.fN; i++) 
254     if (fActive[i] == vol) return kTRUE;
255   return kFALSE;
256 }
257   
258 //____________________________________________________________________
259 AliFMDDetector*
260 AliFMDGeometry::GetDetector(Int_t i) const
261 {
262   // Get the ith detector.   i should be one of 1, 2, or 3.  If an
263   // invalid value is passed, 0 (NULL) is returned. 
264   switch (i) {
265   case 1: return fUseFMD1 ? static_cast<AliFMDDetector*>(fFMD1) : 0;
266   case 2: return fUseFMD2 ? static_cast<AliFMDDetector*>(fFMD2) : 0;
267   case 3: return fUseFMD3 ? static_cast<AliFMDDetector*>(fFMD3) : 0;
268   }
269   return 0;
270 }
271
272 //____________________________________________________________________
273 AliFMDRing*
274 AliFMDGeometry::GetRing(Char_t i) const
275 {
276   // Get the ith ring.  i should be one of 'I' or 'O' (case
277   // insensitive).  If an invalid parameter is passed, 0 (NULL) is
278   // returned. 
279   switch (i) {
280   case 'I':
281   case 'i': return fInner;
282   case 'O':
283   case 'o': return fOuter;
284   }
285   return 0;
286 }
287
288 //____________________________________________________________________
289 void
290 AliFMDGeometry::Enable(Int_t i)
291 {
292   // Enable the ith detector.  i should be one of 1, 2, or 3
293   switch (i) {
294   case 1: fUseFMD1 = kTRUE; break;
295   case 2: fUseFMD2 = kTRUE; break;
296   case 3: fUseFMD3 = kTRUE; break;
297   }
298 }
299
300 //____________________________________________________________________
301 void
302 AliFMDGeometry::Disable(Int_t i)
303 {
304   // Disable the ith detector.  i should be one of 1, 2, or 3
305   switch (i) {
306   case 1: fUseFMD1 = kFALSE; break;
307   case 2: fUseFMD2 = kFALSE; break;
308   case 3: fUseFMD3 = kFALSE; break;
309   }
310 }
311
312 //____________________________________________________________________
313 void
314 AliFMDGeometry::Detector2XYZ(UShort_t  detector, 
315                              Char_t    ring, 
316                              UShort_t  sector, 
317                              UShort_t  strip, 
318                              Double_t& x, 
319                              Double_t& y, 
320                              Double_t& z) const
321 {
322   // Translate detector coordinates (detector, ring, sector, strip) to
323   // spatial coordinates (x, y, z) in the master reference frame of
324   // ALICE. 
325   AliFMDDetector* det = GetDetector(detector);
326   if (!det) { 
327     AliWarning(Form("Unknown detector %d", detector));
328     return;
329   }
330   det->Detector2XYZ(ring, sector, strip, x, y, z);
331 }
332
333 //____________________________________________________________________
334 Bool_t
335 AliFMDGeometry::XYZ2Detector(Double_t  x, 
336                              Double_t  y, 
337                              Double_t  z,
338                              UShort_t& detector, 
339                              Char_t&   ring, 
340                              UShort_t& sector, 
341                              UShort_t& strip) const
342 {
343   // Translate spatial coordinates (x,y,z) in the master reference frame of
344   // ALICE to the detector coordinates (detector, ring, sector,
345   // strip).  Note, that if this method is to be used in
346   // reconstruction or the like, then the input z-coordinate should be
347   // corrected for the events interactions points z-coordinate, like 
348   // geom->XYZ2Detector(x,y,z-ipz,d,r,s,t);
349   AliFMDDetector* det = 0;
350   detector = 0;
351   for (int i = 1; i <= 3; i++) {
352     det = GetDetector(i);
353     if (!det) continue;
354     if (det->XYZ2Detector(x, y, z, ring, sector, strip)) {
355       detector = det->GetId();
356       return kTRUE;
357     }
358   }
359   return kFALSE;
360 }
361
362
363 //____________________________________________________________________
364 void
365 AliFMDGeometry::GetGlobal(const AliRecPoint* p, 
366                           TVector3& pos, 
367                           TMatrixF& /* mat */) const 
368 {
369   // Get the global coordinates cooresponding to the reconstructed
370   // point p.  The coordiates is returned in the 3-vector pos passed
371   // to this member function.  The matrix mat is used for rotations. 
372   GetGlobal(p, pos);
373 }
374
375 //____________________________________________________________________
376 void
377 AliFMDGeometry::GetGlobal(const AliRecPoint* p, TVector3& pos) const 
378 {
379   // Get the global coordinates cooresponding to the reconstructed
380   // point p.  The coordiates is returned in the 3-vector pos passed
381   // to this member function. Note, as AliRecPoint only has places for
382   // 3 indicies, it is assumed that the ring hit is an inner ring -
383   // which obviously needn't be the case. This makes the member
384   // function pretty darn useless. 
385   // FIXME: Implement this function to work with outer rings too. 
386   Double_t x, y, z;
387   TVector3 local;
388   p->GetLocalPosition(local);
389   UShort_t detector = UShort_t(local.X());
390   UShort_t sector   = UShort_t(local.Y());
391   UShort_t strip    = UShort_t(local.Z());
392   Detector2XYZ(detector, 'I', sector, strip, x, y, z);
393   pos.SetXYZ(x, y, z);
394 }
395
396 //____________________________________________________________________
397 Bool_t
398 AliFMDGeometry::Impact(const TParticle* /* particle */) const 
399
400   // Return true, if the particle will hit the active detector
401   // elements, and false if not.  Should be used for fast
402   // simulations.  Note, that the function currently return false
403   // always.  
404   // FIXME: Implement this function. 
405   return kFALSE; 
406 }
407
408 //____________________________________________________________________  
409 void  
410 AliFMDGeometry::SetAlignableVolumes() const
411 {
412   // Declare alignable volumes
413   for (Int_t d = 1; d <= 3; d++) 
414     if (GetDetector(d)) GetDetector(d)->SetAlignableVolumes();
415 }
416
417
418 //____________________________________________________________________  
419 void  
420 AliFMDGeometry::ExtractGeomInfo()
421 {
422   // Check the volume depth of some nodes, get the active volume
423   // numbers, and so forth. 
424   // 
425   // TODO: Here, we should actually also get the parameters of the
426   // shapes, like the verticies of the polygon shape that makes up the
427   // silicon sensor, the strip pitch, the ring radii, the z-positions,
428   // and so on - that is, all the geometric information we need for
429   // futher processing, such as simulation, digitization,
430   // reconstruction, etc. 
431   Int_t detectorDepth = FindNodeDepth("F1MT_1", "ALIC");
432   Int_t ringDepth     = FindNodeDepth(Form("FITV_%d", int('I')), "ALIC");
433   Int_t moduleDepth   = FindNodeDepth("FIBH_0", "ALIC");
434   Int_t sectorDepth   = FindNodeDepth("FISC_1", "ALIC");
435   fActive.Set(0);
436   fActive.Reset(-1);
437   AliFMDDebug(1, ("Geometry depths:\n"
438                    "   Sector:     %d\n"
439                    "   Module:     %d\n"
440                    "   Ring:       %d\n"
441                    "   Detector:   %d", 
442                    sectorDepth, moduleDepth, ringDepth, detectorDepth));
443   if (sectorDepth < 0 && moduleDepth < 0) {
444     fDetailed    = kFALSE;
445     fSectorOff   = -1;
446     fModuleOff   = -1;
447     fRingOff     = 0;
448     fDetectorOff = (ringDepth - detectorDepth);
449     TGeoVolume* actiVol = gGeoManager->GetVolume("FIAC");
450     TGeoVolume* actoVol = gGeoManager->GetVolume("FOAC");
451     if (actiVol) AddActive(actiVol->GetNumber());
452     if (actiVol) AddActive(actoVol->GetNumber());
453   }
454   else if (sectorDepth < 0) {
455     fDetailed    = kFALSE;
456     fSectorOff   = -1;
457     fModuleOff   = 1;
458     fRingOff     = (moduleDepth - ringDepth) + 1;
459     fDetectorOff = (moduleDepth - detectorDepth) + 1;
460     TGeoVolume* modiVol = gGeoManager->GetVolume("FIMO");
461     TGeoVolume* modoVol = gGeoManager->GetVolume("FOMO");
462     if (modiVol) AddActive(modiVol->GetNumber());
463     if (modoVol) AddActive(modoVol->GetNumber());
464   }
465   else {
466     Int_t stripDepth    = FindNodeDepth("FIST_1", "ALIC");
467     fDetailed    = kTRUE;
468     fSectorOff   = (stripDepth - sectorDepth);
469     fModuleOff   = (moduleDepth >= 0 ? (stripDepth - moduleDepth) : -1);
470     fRingOff     = (stripDepth - ringDepth);
471     fDetectorOff = (stripDepth - detectorDepth );
472     TGeoVolume* striVol = gGeoManager->GetVolume("FIST");
473     TGeoVolume* stroVol = gGeoManager->GetVolume("FOST");
474     if (striVol) AddActive(striVol->GetNumber());
475     if (stroVol) AddActive(stroVol->GetNumber());
476   }    
477   AliFMDDebug(1, ("Geometry offsets:\n"
478                    "   Sector:     %d\n"
479                    "   Module:     %d\n"
480                    "   Ring:       %d\n"
481                    "   Detector:   %d", 
482                    fSectorOff, fModuleOff, fRingOff, fDetectorOff));
483 }
484
485   
486 //____________________________________________________________________  
487 static Int_t 
488 CheckNodes(TGeoNode* node, const char* name, Int_t& lvl)
489 {
490   // If there's no node here. 
491   if (!node) return -1;
492   // Check if it this one 
493   TString sname(name);
494   if (sname == node->GetName()) return lvl;
495
496   // Check if the node is an immediate daugther 
497   TObjArray* nodes = node->GetNodes();
498   if (!nodes) return -1;
499   // Increase the level, and search immediate sub nodes. 
500   lvl++;
501   TGeoNode*  found = static_cast<TGeoNode*>(nodes->FindObject(name));
502   if (found) return lvl;
503
504   // Check the sub node, if any of their sub-nodes match.
505   for (Int_t i = 0; i < nodes->GetEntries(); i++) {
506     TGeoNode* sub = static_cast<TGeoNode*>(nodes->At(i));
507     if (!sub) continue;
508     // Recurive check 
509     if (CheckNodes(sub, name, lvl) >= 0) return lvl;
510   }
511   // If not found, decrease the level 
512   lvl--;
513   return -1;
514 }
515 //____________________________________________________________________  
516 Int_t 
517 FindNodeDepth(const char* name, const char* volname) 
518 {
519   // Find the depth of a node 
520   TGeoVolume* vol  = gGeoManager->GetVolume(volname);
521   if (!vol) {
522     std::cerr << "No top volume defined" << std::endl;
523     return -1;
524   }
525
526   TGeoIterator next(vol);
527   TGeoNode*    node = 0;
528   TString      sName(name);
529   while ((node = next())) { 
530     // std::cout << "Checking node " << node->GetName() << std::endl;
531     if (sName == node->GetName()) { 
532       std::cout << "Found node " << node->GetName() << " at level " 
533                 << next.GetLevel() << std::endl;
534       return next.GetLevel();
535     }
536   }
537   return -1;
538 #if 0
539   TObjArray* nodes = vol->GetNodes();
540   if (!nodes) { 
541     std::cerr << "No nodes in top volume" << std::endl;
542     return -1;
543   }
544   TIter next(nodes);
545   TGeoNode* node = 0;
546   Int_t lvl = 0;
547   while ((node = static_cast<TGeoNode*>(next()))) 
548     if (CheckNodes(node, name, lvl) >= 0) return lvl;
549   return -1;
550 #endif
551 }
552
553 //____________________________________________________________________
554 //
555 // EOF
556 //