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 //