[u/mrichter/AliRoot.git] / FMD / AliFMDGeometry.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
/* $Id$ */
/** @file    AliFMDGeometry.cxx
    @author  Christian Holm Christensen <cholm@nbi.dk>
    @date    Mon Mar 27 12:40:37 2006
    @brief   Geometry mananger for the FMD
*/
//____________________________________________________________________
//                                                                          
// Forward Multiplicity Detector based on Silicon wafers. 
//
// This class is a singleton that handles the geometry parameters of
// the FMD detectors.  
//                                                       
// The actual code is done by various separate classes.   Below is
// diagram showing the relationship between the various FMD classes
// that handles the geometry 
//
//                               +------------+ 
//                            +- | AliFMDRing |
// 			   2  |  +------------+
//      +----------------+<>--+        |				
//      | AliFMDGeometry |             ^                       	
//      +----------------+<>--+        V 1..2                     	
//           		   3  | +----------------+ 		
//            		      +-| AliFMDDetector | 		
//             		        +----------------+		
//                                     ^
//                                     |
//                       +-------------+-------------+
//                       |             |             |	      
//                  +---------+   +---------+   +---------+
//                  | AliFMD1 |   | AliFMD2 |   | AliFMD3 |
//                  +---------+   +---------+   +---------+
//      
//
// *  AliFMDRing 
//    This class contains all stuff needed to do with a ring.  It's
//    used by the AliFMDDetector objects to instantise inner and
//    outer rings.  The AliFMDRing objects are shared by the
//    AliFMDDetector objects, and owned by the AliFMDv1 object. 
//
// *  AliFMD1, AliFMD2, and AliFMD3 
//    These are specialisation of AliFMDDetector, that contains the
//    particularities of each of the sub-detector system.  It is
//    envisioned that the classes should also define the support
//    volumes and material for each of the detectors.                          
//                                                                          
//
#include "AliFMDGeometry.h"	// ALIFMDGEOMETRY_H
#include "AliFMDRing.h"		// ALIFMDRING_H
#include "AliFMD1.h"		// ALIFMD1_H
#include "AliFMD2.h"		// ALIFMD2_H
#include "AliFMD3.h"		// ALIFMD2_H
#include "AliRecPoint.h"	// ALIRECPOINT_H
#include "AliFMDDebug.h"		   // ALILOG_H
#include <TVector3.h>           // ROOT_TVector3
// #include <TMatrix.h>            // ROOT_TMatrix
// #include <TParticle.h>          // ROOT_TParticle
#include <Riostream.h>
#include "AliFMDGeometryBuilder.h"
// #include <TArrayI.h>
#include <TGeoManager.h>
#include <TGeoVolume.h>
#include <TGeoNode.h>
#include <TMath.h>
static Int_t FindNodeDepth(const char* name, const char* volname);


//====================================================================
ClassImp(AliFMDGeometry)
#if 0
  ; // This is here to keep Emacs for indenting the next line
#endif

//____________________________________________________________________
AliFMDGeometry* AliFMDGeometry::fgInstance = 0;

//____________________________________________________________________
AliFMDGeometry* 
AliFMDGeometry::Instance() 
{
  // Return (newly created) singleton instance 
  if (!fgInstance) fgInstance = new AliFMDGeometry;
  return fgInstance;
}

//____________________________________________________________________
AliFMDGeometry::AliFMDGeometry() 
  : AliGeometry("FMD", "Forward multiplicity"), 
    fIsInitialized(kFALSE), 
    fInner(0),
    fOuter(0),
    fFMD1(0),
    fFMD2(0),
    fFMD3(0),
    fUseFMD1(kFALSE),
    fUseFMD2(kFALSE),
    fUseFMD3(kFALSE),
    fIsInitTrans(kFALSE),
    fBuilder(0),
    fDetectorOff(0),
    fModuleOff(0),  
    fRingOff(0),
    fSectorOff(0),
    fActive(2),
    fDetailed(kFALSE),       
    fUseAssembly(kFALSE)
{
  // PROTECTED
  // Default constructor 
  fUseFMD1     = kTRUE;
  fUseFMD2     = kTRUE;
  fUseFMD3     = kTRUE;  
  fDetailed    = kTRUE;
  fUseAssembly = kTRUE;
  fInner = new AliFMDRing('I');
  fOuter = new AliFMDRing('O');
  fFMD1  = new AliFMD1(fInner);
  fFMD2  = new AliFMD2(fInner, fOuter);
  fFMD3  = new AliFMD3(fInner, fOuter);
  fIsInitialized = kFALSE;
  fActive.Reset(-1);
}

//____________________________________________________________________
AliFMDGeometry::AliFMDGeometry(const AliFMDGeometry& other) 
  : AliGeometry(other),
    fIsInitialized(other.fIsInitialized),
    fInner(other.fInner), 
    fOuter(other.fOuter), 
    fFMD1(other.fFMD1), 
    fFMD2(other.fFMD2), 
    fFMD3(other.fFMD3), 
    fUseFMD1(other.fUseFMD1), 
    fUseFMD2(other.fUseFMD2), 
    fUseFMD3(other.fUseFMD3), 
    fIsInitTrans(other.fIsInitTrans),
    fBuilder(other.fBuilder),
    fDetectorOff(other.fDetectorOff),
    fModuleOff(other.fModuleOff),  
    fRingOff(other.fRingOff),
    fSectorOff(other.fSectorOff),
    fActive(other.fActive),
    fDetailed(other.fDetailed),
    fUseAssembly(other.fUseAssembly)
{
  // PROTECTED
  // Copy constructor
}


//____________________________________________________________________
AliFMDGeometry&
AliFMDGeometry::operator=(const AliFMDGeometry& other) 
{
  // PROTECTED
  // Assignment operator 
  fUseFMD1		= other.fUseFMD1; 
  fUseFMD2		= other.fUseFMD2; 
  fUseFMD3		= other.fUseFMD3; 
  fFMD1			= other.fFMD1; 
  fFMD2			= other.fFMD2; 
  fFMD3			= other.fFMD3; 
  fInner		= other.fInner; 
  fOuter		= other.fOuter; 
  fIsInitialized	= other.fIsInitialized;
  return *this;
}

//____________________________________________________________________
void
AliFMDGeometry::Init()
{
  // Initialize the the singleton if not done so already 
  if (fIsInitialized) return;
  fInner->Init();
  fOuter->Init();
  fFMD1->Init();
  fFMD2->Init();
  fFMD3->Init();
}

//____________________________________________________________________
void
AliFMDGeometry::InitTransformations(Bool_t force)
{
  // Find all local <-> global transforms 
  if (force) fIsInitTrans = kFALSE;
  if (fIsInitTrans) return; 
  if (!gGeoManager) {
    AliError("No TGeoManager defined");
    return;
  }
  AliFMDDebug(1, ("Initialising transforms for FMD geometry"));
  if (fFMD1) fFMD1->InitTransformations();
  if (fFMD2) fFMD2->InitTransformations();
  if (fFMD3) fFMD3->InitTransformations();
  fIsInitTrans = kTRUE;
}

//____________________________________________________________________
void
AliFMDGeometry::Build()
{
  // Build the geometry 
  if (!fBuilder) fBuilder = new AliFMDGeometryBuilder(fDetailed);
  fBuilder->SetDetailed(fDetailed);
  fBuilder->UseAssembly(fUseAssembly);
  fBuilder->Exec();
}

//____________________________________________________________________
void
AliFMDGeometry::SetActive(Int_t* active, Int_t n) 
{
  // Set active volumes 
  fActive.Set(n);
  for (Int_t i = 0; i < n; i++) { 
    AliFMDDebug(1, ("Active vol id # %d: %d", i, active[i]));
    fActive[i] = active[i];
  }
}

//____________________________________________________________________
void
AliFMDGeometry::AddActive(Int_t active)
{
  // Add an active volume 
  Int_t n = fActive.fN;
  fActive.Set(n+1);
  fActive[n] = active;
}

//____________________________________________________________________
Bool_t
AliFMDGeometry::IsActive(Int_t vol) const
{
  // Check if a volume is active 
  for (Int_t i = 0; i < fActive.fN; i++) 
    if (fActive[i] == vol) return kTRUE;
  return kFALSE;
}
  
//____________________________________________________________________
AliFMDDetector*
AliFMDGeometry::GetDetector(Int_t i) const
{
  // Get the ith detector.   i should be one of 1, 2, or 3.  If an
  // invalid value is passed, 0 (NULL) is returned. 
  switch (i) {
  case 1: return fUseFMD1 ? static_cast<AliFMDDetector*>(fFMD1) : 0;
  case 2: return fUseFMD2 ? static_cast<AliFMDDetector*>(fFMD2) : 0;
  case 3: return fUseFMD3 ? static_cast<AliFMDDetector*>(fFMD3) : 0;
  }
  return 0;
}

//____________________________________________________________________
AliFMDRing*
AliFMDGeometry::GetRing(Char_t i) const
{
  // Get the ith ring.  i should be one of 'I' or 'O' (case
  // insensitive).  If an invalid parameter is passed, 0 (NULL) is
  // returned. 
  switch (i) {
  case 'I':
  case 'i': return fInner;
  case 'O':
  case 'o': return fOuter;
  }
  return 0;
}

//____________________________________________________________________
void
AliFMDGeometry::Enable(Int_t i)
{
  // Enable the ith detector.  i should be one of 1, 2, or 3
  switch (i) {
  case 1: fUseFMD1 = kTRUE; break;
  case 2: fUseFMD2 = kTRUE; break;
  case 3: fUseFMD3 = kTRUE; break;
  }
}

//____________________________________________________________________
void
AliFMDGeometry::Disable(Int_t i)
{
  // Disable the ith detector.  i should be one of 1, 2, or 3
  switch (i) {
  case 1: fUseFMD1 = kFALSE; break;
  case 2: fUseFMD2 = kFALSE; break;
  case 3: fUseFMD3 = kFALSE; break;
  }
}

//____________________________________________________________________
void
AliFMDGeometry::Detector2XYZ(UShort_t  detector, 
			     Char_t    ring, 
			     UShort_t  sector, 
			     UShort_t  strip, 
			     Double_t& x, 
			     Double_t& y, 
			     Double_t& z) const
{
  // Translate detector coordinates (detector, ring, sector, strip) to
  // spatial coordinates (x, y, z) in the master reference frame of
  // ALICE. 
  AliFMDDetector* det = GetDetector(detector);
  if (!det) { 
    AliWarning(Form("Unknown detector %d", detector));
    return;
  }
  det->Detector2XYZ(ring, sector, strip, x, y, z);
}

//____________________________________________________________________
Bool_t
AliFMDGeometry::XYZ2Detector(Double_t  x, 
			     Double_t  y, 
			     Double_t  z,
			     UShort_t& detector, 
			     Char_t&   ring, 
			     UShort_t& sector, 
			     UShort_t& strip) const
{
  // Translate spatial coordinates (x,y,z) in the master reference frame of
  // ALICE to the detector coordinates (detector, ring, sector,
  // strip).  Note, that if this method is to be used in
  // reconstruction or the like, then the input z-coordinate should be
  // corrected for the events interactions points z-coordinate, like 
  // geom->XYZ2Detector(x,y,z-ipz,d,r,s,t);
  AliFMDDetector* det = 0;
  detector = 0;
  for (int i = 1; i <= 3; i++) {
    det = GetDetector(i);
    if (!det) continue;
    if (det->XYZ2Detector(x, y, z, ring, sector, strip)) {
      detector = det->GetId();
      return kTRUE;
    }
  }
  return kFALSE;
}

//____________________________________________________________________
Bool_t
AliFMDGeometry::XYZ2REtaPhiTheta(Double_t  x,   Double_t y, 
				 Double_t  z, 
				 Double_t& r,   Double_t& eta, 
				 Double_t& phi, Double_t& theta)
{
  
  // Service function to convert Cartisean XYZ to r, eta, phi, and theta.   
  // 
  // Note, that the z input should be corrected for the vertex location 
  // if needed.
  //
  //     x      Cartisean X coordinate
  //     y      Cartisean Y coordinate 
  //     z      Cartisean Z coordinate 
  //     r      On return, the radius
  //     eta    On return, the pseudo-rapidity
  //     phi    On return, the azimuthal angle
  //     theta  On return, the polar angle;
  // 
  // Return:
  //     kFALSE in case of problems. 

  if (x == 0 && y == 0 && z == 0) return kFALSE;
  
  // Correct for vertex offset. 
  phi   =  TMath::ATan2(y, x);
  r     =  TMath::Sqrt(y * y + x * x);
  theta =  TMath::ATan2(r, z);
  eta   = -TMath::Log(TMath::Tan(theta / 2));

  return kTRUE;
}


//____________________________________________________________________
void
AliFMDGeometry::GetGlobal(const AliRecPoint* p, 
			  TVector3& pos, 
			  TMatrixF& /* mat */) const 
{
  // Get the global coordinates cooresponding to the reconstructed
  // point p.  The coordiates is returned in the 3-vector pos passed
  // to this member function.  The matrix mat is used for rotations. 
  GetGlobal(p, pos);
}

//____________________________________________________________________
void
AliFMDGeometry::GetGlobal(const AliRecPoint* p, TVector3& pos) const 
{
  // Get the global coordinates cooresponding to the reconstructed
  // point p.  The coordiates is returned in the 3-vector pos passed
  // to this member function. Note, as AliRecPoint only has places for
  // 3 indicies, it is assumed that the ring hit is an inner ring -
  // which obviously needn't be the case. This makes the member
  // function pretty darn useless. 
  // FIXME: Implement this function to work with outer rings too. 
  Double_t x, y, z;
  TVector3 local;
  p->GetLocalPosition(local);
  UShort_t detector = UShort_t(local.X());
  UShort_t sector   = UShort_t(local.Y());
  UShort_t strip    = UShort_t(local.Z());
  Detector2XYZ(detector, 'I', sector, strip, x, y, z);
  pos.SetXYZ(x, y, z);
}

//____________________________________________________________________
Bool_t
AliFMDGeometry::Impact(const TParticle* /* particle */) const 
{ 
  // Return true, if the particle will hit the active detector
  // elements, and false if not.  Should be used for fast
  // simulations.  Note, that the function currently return false
  // always.  
  // FIXME: Implement this function. 
  return kFALSE; 
}

//____________________________________________________________________	
void  
AliFMDGeometry::SetAlignableVolumes() const
{
  // Declare alignable volumes
  for (Int_t d = 1; d <= 3; d++) 
    if (GetDetector(d)) GetDetector(d)->SetAlignableVolumes();
}


//____________________________________________________________________	
void  
AliFMDGeometry::ExtractGeomInfo()
{
  // Check the volume depth of some nodes, get the active volume
  // numbers, and so forth. 
  // 
  // TODO: Here, we should actually also get the parameters of the
  // shapes, like the verticies of the polygon shape that makes up the
  // silicon sensor, the strip pitch, the ring radii, the z-positions,
  // and so on - that is, all the geometric information we need for
  // futher processing, such as simulation, digitization,
  // reconstruction, etc. 
  Int_t detectorDepth = FindNodeDepth("F1MT_1", "ALIC");
  Int_t ringDepth     = FindNodeDepth(Form("FITV_%d", int('I')), "ALIC");
  Int_t moduleDepth   = FindNodeDepth("FIBH_0", "ALIC");
  Int_t sectorDepth   = FindNodeDepth("FISC_1", "ALIC");
  fActive.Set(0);
  fActive.Reset(-1);
  AliFMDDebug(1, ("Geometry depths:\n"
		   "   Sector:     %d\n"
		   "   Module:     %d\n"
		   "   Ring:       %d\n"
		   "   Detector:   %d", 
		   sectorDepth, moduleDepth, ringDepth, detectorDepth));
  if (sectorDepth < 0 && moduleDepth < 0) {
    fDetailed    = kFALSE;
    fSectorOff   = -1;
    fModuleOff   = -1;
    fRingOff     = 0;
    fDetectorOff = (ringDepth - detectorDepth);
    TGeoVolume* actiVol = gGeoManager->GetVolume("FIAC");
    TGeoVolume* actoVol = gGeoManager->GetVolume("FOAC");
    if (actiVol) AddActive(actiVol->GetNumber());
    if (actiVol) AddActive(actoVol->GetNumber());
  }
  else if (sectorDepth < 0) {
    fDetailed    = kFALSE;
    fSectorOff   = -1;
    fModuleOff   = 1;
    fRingOff     = (moduleDepth - ringDepth) + 1;
    fDetectorOff = (moduleDepth - detectorDepth) + 1;
    TGeoVolume* modiVol = gGeoManager->GetVolume("FIMO");
    TGeoVolume* modoVol = gGeoManager->GetVolume("FOMO");
    if (modiVol) AddActive(modiVol->GetNumber());
    if (modoVol) AddActive(modoVol->GetNumber());
  }
  else {
    Int_t stripDepth    = FindNodeDepth("FIST_1", "ALIC");
    fDetailed    = kTRUE;
    fSectorOff   = (stripDepth - sectorDepth);
    fModuleOff   = (moduleDepth >= 0 ? (stripDepth - moduleDepth) : -1);
    fRingOff     = (stripDepth - ringDepth);
    fDetectorOff = (stripDepth - detectorDepth );
    TGeoVolume* striVol = gGeoManager->GetVolume("FIST");
    TGeoVolume* stroVol = gGeoManager->GetVolume("FOST");
    if (striVol) AddActive(striVol->GetNumber());
    if (stroVol) AddActive(stroVol->GetNumber());
  }    
  AliFMDDebug(1, ("Geometry offsets:\n"
		   "   Sector:     %d\n"
		   "   Module:     %d\n"
		   "   Ring:       %d\n"
		   "   Detector:   %d", 
		   fSectorOff, fModuleOff, fRingOff, fDetectorOff));
}

  
//____________________________________________________________________	
static Int_t 
CheckNodes(TGeoNode* node, const char* name, Int_t& lvl)
{
  // If there's no node here. 
  if (!node) return -1;
  // Check if it this one 
  TString sname(name);
  if (sname == node->GetName()) return lvl;

  // Check if the node is an immediate daugther 
  TObjArray* nodes = node->GetNodes();
  if (!nodes) return -1;
  // Increase the level, and search immediate sub nodes. 
  lvl++;
  TGeoNode*  found = static_cast<TGeoNode*>(nodes->FindObject(name));
  if (found) return lvl;

  // Check the sub node, if any of their sub-nodes match.
  for (Int_t i = 0; i < nodes->GetEntries(); i++) {
    TGeoNode* sub = static_cast<TGeoNode*>(nodes->At(i));
    if (!sub) continue;
    // Recurive check 
    if (CheckNodes(sub, name, lvl) >= 0) return lvl;
  }
  // If not found, decrease the level 
  lvl--;
  return -1;
}
//____________________________________________________________________	
Int_t 
FindNodeDepth(const char* name, const char* volname) 
{
  // Find the depth of a node 
  TGeoVolume* vol  = gGeoManager->GetVolume(volname);
  if (!vol) {
    std::cerr << "No top volume defined" << std::endl;
    return -1;
  }

  TGeoIterator next(vol);
  TGeoNode*    node = 0;
  TString      sName(name);
  while ((node = next())) { 
    if (sName == node->GetName()) { 
      //std::cout << "Found node " << node->GetName() << " at level " 
      //		<< next.GetLevel() << std::endl;
      return next.GetLevel();
    }
  }
  return -1;
#if 0
  TObjArray* nodes = vol->GetNodes();
  if (!nodes) { 
    std::cerr << "No nodes in top volume" << std::endl;
    return -1;
  }
  TIter next(nodes);
  TGeoNode* node = 0;
  Int_t lvl = 0;
  while ((node = static_cast<TGeoNode*>(next()))) 
    if (CheckNodes(node, name, lvl) >= 0) return lvl;
  return -1;
#endif
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
1a1fdef7	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	**************************************************************************/
1a1fdef7	15	/* $Id$ */
c2fc1258	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	*/
1a1fdef7	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
f95a63c4	69	#include "AliFMDDebug.h" // ALILOG_H
1a1fdef7	70	#include <TVector3.h> // ROOT_TVector3
02a27b50	71	// #include <TMatrix.h> // ROOT_TMatrix
02a27b50	72	// #include <TParticle.h> // ROOT_TParticle
1a1fdef7	73	#include <Riostream.h>
54e415a8	74	#include "AliFMDGeometryBuilder.h"
02a27b50	75	// #include <TArrayI.h>
54e415a8	76	#include <TGeoManager.h>
	77	#include <TGeoVolume.h>
	78	#include <TGeoNode.h>
75609cab	79	#include <TMath.h>
54e415a8	80	static Int_t FindNodeDepth(const char* name, const char* volname);
54e415a8	81
1a1fdef7	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	{
69b696b9	96	// Return (newly created) singleton instance
1a1fdef7	97	if (!fgInstance) fgInstance = new AliFMDGeometry;
	98	return fgInstance;
	99	}
	100
	101	//____________________________________________________________________
	102	AliFMDGeometry::AliFMDGeometry()
a1f80595	103	: AliGeometry("FMD", "Forward multiplicity"),
b5ee4425	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),
352482dc	113	fIsInitTrans(kFALSE),
b5ee4425	114	fBuilder(0),
	115	fDetectorOff(0),
	116	fModuleOff(0),
	117	fRingOff(0),
	118	fSectorOff(0),
	119	fActive(2),
	120	fDetailed(kFALSE),
	121	fUseAssembly(kFALSE)
1a1fdef7	122	{
69b696b9	123	// PROTECTED
69b696b9	124	// Default constructor
54e415a8	125	fUseFMD1 = kTRUE;
	126	fUseFMD2 = kTRUE;
	127	fUseFMD3 = kTRUE;
54e415a8	128	fDetailed = kTRUE;
54e415a8	129	fUseAssembly = kTRUE;
1a1fdef7	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;
54e415a8	136	fActive.Reset(-1);
1a1fdef7	137	}
1a1fdef7	138
69b696b9	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),
bf000c32	150	fUseFMD3(other.fUseFMD3),
352482dc	151	fIsInitTrans(other.fIsInitTrans),
b5ee4425	152	fBuilder(other.fBuilder),
	153	fDetectorOff(other.fDetectorOff),
	154	fModuleOff(other.fModuleOff),
	155	fRingOff(other.fRingOff),
	156	fSectorOff(other.fSectorOff),
bf000c32	157	fActive(other.fActive),
54e415a8	158	fDetailed(other.fDetailed),
54e415a8	159	fUseAssembly(other.fUseAssembly)
69b696b9	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
1a1fdef7	185	//____________________________________________________________________
	186	void
	187	AliFMDGeometry::Init()
	188	{
69b696b9	189	// Initialize the the singleton if not done so already
1a1fdef7	190	if (fIsInitialized) return;
	191	fInner->Init();
	192	fOuter->Init();
	193	fFMD1->Init();
	194	fFMD2->Init();
	195	fFMD3->Init();
	196	}
	197
bf000c32	198	//____________________________________________________________________
bf000c32	199	void
352482dc	200	AliFMDGeometry::InitTransformations(Bool_t force)
bf000c32	201	{
02a27b50	202	// Find all local <-> global transforms
352482dc	203	if (force) fIsInitTrans = kFALSE;
352482dc	204	if (fIsInitTrans) return;
bf000c32	205	if (!gGeoManager) {
	206	AliError("No TGeoManager defined");
	207	return;
	208	}
352482dc	209	AliFMDDebug(1, ("Initialising transforms for FMD geometry"));
bf000c32	210	if (fFMD1) fFMD1->InitTransformations();
	211	if (fFMD2) fFMD2->InitTransformations();
	212	if (fFMD3) fFMD3->InitTransformations();
352482dc	213	fIsInitTrans = kTRUE;
bf000c32	214	}
bf000c32	215
54e415a8	216	//____________________________________________________________________
	217	void
	218	AliFMDGeometry::Build()
	219	{
02a27b50	220	// Build the geometry
54e415a8	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	{
02a27b50	231	// Set active volumes
54e415a8	232	fActive.Set(n);
088f8e79	233	for (Int_t i = 0; i < n; i++) {
f95a63c4	234	AliFMDDebug(1, ("Active vol id # %d: %d", i, active[i]));
088f8e79	235	fActive[i] = active[i];
088f8e79	236	}
54e415a8	237	}
	238
	239	//____________________________________________________________________
	240	void
	241	AliFMDGeometry::AddActive(Int_t active)
	242	{
02a27b50	243	// Add an active volume
54e415a8	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	{
02a27b50	253	// Check if a volume is active
54e415a8	254	for (Int_t i = 0; i < fActive.fN; i++)
	255	if (fActive[i] == vol) return kTRUE;
	256	return kFALSE;
	257	}
54e415a8	258
1a1fdef7	259	//____________________________________________________________________
	260	AliFMDDetector*
	261	AliFMDGeometry::GetDetector(Int_t i) const
	262	{
69b696b9	263	// Get the ith detector. i should be one of 1, 2, or 3. If an
69b696b9	264	// invalid value is passed, 0 (NULL) is returned.
1a1fdef7	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	{
69b696b9	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.
1a1fdef7	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	{
69b696b9	293	// Enable the ith detector. i should be one of 1, 2, or 3
1a1fdef7	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	{
69b696b9	305	// Disable the ith detector. i should be one of 1, 2, or 3
1a1fdef7	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,
69b696b9	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
1a1fdef7	322	{
69b696b9	323	// Translate detector coordinates (detector, ring, sector, strip) to
	324	// spatial coordinates (x, y, z) in the master reference frame of
	325	// ALICE.
1a1fdef7	326	AliFMDDetector* det = GetDetector(detector);
bf000c32	327	if (!det) {
	328	AliWarning(Form("Unknown detector %d", detector));
	329	return;
	330	}
1a1fdef7	331	det->Detector2XYZ(ring, sector, strip, x, y, z);
	332	}
	333
54e415a8	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
75609cab	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
1a1fdef7	398
	399	//____________________________________________________________________
	400	void
	401	AliFMDGeometry::GetGlobal(const AliRecPoint* p,
	402	TVector3& pos,
e8d02863	403	TMatrixF& /* mat */) const
1a1fdef7	404	{
69b696b9	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.
1a1fdef7	408	GetGlobal(p, pos);
	409	}
	410
	411	//____________________________________________________________________
	412	void
	413	AliFMDGeometry::GetGlobal(const AliRecPoint* p, TVector3& pos) const
	414	{
69b696b9	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.
1a1fdef7	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	{
69b696b9	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.
1a1fdef7	441	return kFALSE;
	442	}
	443
9de78b35	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
54e415a8	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.
53a1201d	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);
54e415a8	472	fActive.Reset(-1);
f95a63c4	473	AliFMDDebug(1, ("Geometry depths:\n"
54e415a8	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	}
f95a63c4	513	AliFMDDebug(1, ("Geometry offsets:\n"
54e415a8	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	{
02a27b50	555	// Find the depth of a node
54e415a8	556	TGeoVolume* vol = gGeoManager->GetVolume(volname);
	557	if (!vol) {
	558	std::cerr << "No top volume defined" << std::endl;
	559	return -1;
	560	}
53a1201d	561
	562	TGeoIterator next(vol);
	563	TGeoNode* node = 0;
	564	TString sName(name);
	565	while ((node = next())) {
53a1201d	566	if (sName == node->GetName()) {
2af4b50b	567	//std::cout << "Found node " << node->GetName() << " at level "
2af4b50b	568	// << next.GetLevel() << std::endl;
53a1201d	569	return next.GetLevel();
	570	}
	571	}
	572	return -1;
	573	#if 0
54e415a8	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;
53a1201d	585	#endif
54e415a8	586	}
54e415a8	587
1a1fdef7	588	//____________________________________________________________________
	589	//
	590	// EOF
	591	//