]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSGeometry.cxx
Added method to calculate absID from position in PHOS local coordinate system
[u/mrichter/AliRoot.git] / PHOS / AliPHOSGeometry.cxx
index 5e0b93df7288756d5452c6533336792c00d17ef4..c1bf473819ffc2ae8e926dcbe6748d21f1357d0e 100644 (file)
 
 #include "TVector3.h"
 #include "TRotation.h" 
-#include "TFolder.h" 
-#include "TROOT.h" 
+#include "TParticle.h"
 
 // --- Standard library ---
 
-#include <iostream.h>
-#include <stdlib.h>
-
 // --- AliRoot header files ---
 
 #include "AliPHOSGeometry.h"
 #include "AliPHOSEMCAGeometry.h" 
 #include "AliPHOSRecPoint.h"
-#include "AliConst.h"
 
 ClassImp(AliPHOSGeometry) ;
 
@@ -51,6 +46,18 @@ ClassImp(AliPHOSGeometry) ;
 AliPHOSGeometry * AliPHOSGeometry::fgGeom = 0 ;
 Bool_t            AliPHOSGeometry::fgInit = kFALSE ;
 
+//____________________________________________________________________________
+AliPHOSGeometry::AliPHOSGeometry() {
+    // default ctor 
+    // must be kept public for root persistency purposes, but should never be called by the outside world
+    fPHOSAngle      = 0 ;
+    fGeometryEMCA   = 0 ;
+    fGeometrySUPP   = 0 ;
+    fGeometryCPV    = 0 ;
+    fgGeom          = 0 ;
+    fRotMatrixArray = 0 ;  
+}  
+
 //____________________________________________________________________________
 AliPHOSGeometry::~AliPHOSGeometry(void)
 {
@@ -66,9 +73,12 @@ void AliPHOSGeometry::Init(void)
 {
   // Initializes the PHOS parameters :
   //  IHEP is the Protvino CPV (cathode pad chambers)
-  //  GPS2 is the Subatech Pre-Shower (two micromegas sandwiching a passive lead converter)
-  //  MIXT 4 PHOS modules withe the IHEP CPV qnd one PHOS module with the Subatche Pre-Shower
   
+  TString test(GetName()) ; 
+  if (test != "IHEP" ) {
+    Fatal("Init", "%s is not a known geometry (choose among IHEP)", test.Data() ) ; 
+  }
+
   fgInit     = kTRUE ; 
   
   fNModules     = 5;
@@ -85,8 +95,8 @@ void AliPHOSGeometry::Init(void)
   Float_t * emcParams = fGeometryEMCA->GetEMCParams() ;
   
   fPHOSParams[0] =  TMath::Max((Double_t)fGeometryCPV->GetCPVBoxSize(0)/2., 
-                              (Double_t)(emcParams[0]*(fGeometryCPV->GetCPVBoxSize(1)+emcParams[3]) - 
-                               emcParams[1]* fGeometryCPV->GetCPVBoxSize(1))/emcParams[3] ) ;
+                              (Double_t)(emcParams[0] - (emcParams[1]-emcParams[0])*
+                                         fGeometryCPV->GetCPVBoxSize(1)/2/emcParams[3]));
   fPHOSParams[1] = emcParams[1] ;
   fPHOSParams[2] = TMath::Max((Double_t)emcParams[2], (Double_t)fGeometryCPV->GetCPVBoxSize(2)/2.);
   fPHOSParams[3] = emcParams[3] + fGeometryCPV->GetCPVBoxSize(1)/2. ;
@@ -102,13 +112,12 @@ void AliPHOSGeometry::Init(void)
   
 }
 
-
 //____________________________________________________________________________
 AliPHOSGeometry *  AliPHOSGeometry::GetInstance() 
 { 
   // Returns the pointer of the unique instance; singleton specific
   
-  return (AliPHOSGeometry *) fgGeom ; 
+  return static_cast<AliPHOSGeometry *>( fgGeom ) ; 
 }
 
 //____________________________________________________________________________
@@ -133,10 +142,8 @@ AliPHOSGeometry *  AliPHOSGeometry::GetInstance(const Text_t* name, const Text_t
     }
   }
   else {
-    if ( strcmp(fgGeom->GetName(), name) != 0 ) {
-      cout << "AliPHOSGeometry <E> : current geometry is " << fgGeom->GetName() << endl
-          << "                      you cannot call     " << name << endl ; 
-    }
+    if ( strcmp(fgGeom->GetName(), name) != 0 ) 
+      ::Error("GetInstance", "Current geometry is %s. You cannot call %s", fgGeom->GetName(), name) ; 
     else
       rv = (AliPHOSGeometry *) fgGeom ; 
   } 
@@ -148,12 +155,11 @@ void AliPHOSGeometry::SetPHOSAngles()
 { 
   // Calculates the position of the PHOS modules in ALICE global coordinate system
   
-  Double_t const kRADDEG = 180.0 / kPI ;
+  Double_t const kRADDEG = 180.0 / TMath::Pi() ;
   Float_t pphi =  2 * TMath::ATan( GetOuterBoxSize(0)  / ( 2.0 * GetIPtoUpperCPVsurface() ) ) ;
   pphi *= kRADDEG ;
   if (pphi > fAngle){ 
-    cout << "AliPHOSGeometry: PHOS modules overlap!\n";
-    cout <<  "pphi = " << pphi << " fAngle  " << fAngle << endl ;
+    Error("SetPHOSAngles", "PHOS modules overlap!\n pphi = %f fAngle = %f", pphi, fAngle);
 
   }
   pphi = fAngle;
@@ -165,7 +171,7 @@ void AliPHOSGeometry::SetPHOSAngles()
 }
 
 //____________________________________________________________________________
-Bool_t AliPHOSGeometry::AbsToRelNumbering(const Int_t AbsId, Int_t * relid) const
+Bool_t AliPHOSGeometry::AbsToRelNumbering(Int_t AbsId, Int_t * relid) const
 {
   // Converts the absolute numbering into the following array/
   //  relid[0] = PHOS Module number 1:fNModules 
@@ -201,7 +207,7 @@ Bool_t AliPHOSGeometry::AbsToRelNumbering(const Int_t AbsId, Int_t * relid) cons
 }
 
 //____________________________________________________________________________  
-void AliPHOSGeometry::EmcModuleCoverage(const Int_t mod, Double_t & tm, Double_t & tM, Double_t & pm, Double_t & pM, Option_t * opt) const 
+void AliPHOSGeometry::EmcModuleCoverage(Int_t mod, Double_t & tm, Double_t & tM, Double_t & pm, Double_t & pM, Option_t * opt) const 
 {
   // calculates the angular coverage in theta and phi of one EMC (=PHOS) module
 
@@ -211,7 +217,7 @@ void AliPHOSGeometry::EmcModuleCoverage(const Int_t mod, Double_t & tm, Double_t
   else if ( opt == Degre() )
     conv = 180. / TMath::Pi() ; 
   else {
-    cout << "<I>  AliPHOSGeometry::EmcXtalCoverage : " << opt << " unknown option; result in radian " << endl ; 
+    Warning("EmcModuleCoverage", "%s unknown option; result in radian", opt) ; 
     conv = 1. ;
       }
 
@@ -246,7 +252,7 @@ void AliPHOSGeometry::EmcXtalCoverage(Double_t & theta, Double_t & phi, Option_t
   else if ( opt == Degre() )
     conv = 180. / TMath::Pi() ; 
   else {
-    cout << "<I>  AliPHOSGeometry::EmcXtalCoverage : " << opt << " unknown option; result in radian " << endl ; 
+    Warning("EmcXtalCoverage", "%s unknown option; result in radian", opt) ;  
     conv = 1. ;
       }
 
@@ -257,7 +263,7 @@ void AliPHOSGeometry::EmcXtalCoverage(Double_t & theta, Double_t & phi, Option_t
  
 
 //____________________________________________________________________________
-void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TMatrix & gmat) const
+void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TMatrix & /*gmat*/) const
 {
   // Calculates the coordinates of a RecPoint and the error matrix in the ALICE global coordinate system
  
@@ -277,7 +283,7 @@ void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TM
     }  
 
   Float_t phi           = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; 
-  Double_t const kRADDEG = 180.0 / kPI ;
+  Double_t const kRADDEG = 180.0 / TMath::Pi() ;
   Float_t rphi          = phi / kRADDEG ; 
   
   TRotation rot ;
@@ -307,7 +313,7 @@ void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos) co
     }  
 
   Float_t phi           = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; 
-  Double_t const kRADDEG = 180.0 / kPI ;
+  Double_t const kRADDEG = 180.0 / TMath::Pi() ;
   Float_t rphi          = phi / kRADDEG ; 
   
   TRotation rot ;
@@ -318,23 +324,24 @@ void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos) co
 }
 
 //____________________________________________________________________________
-void AliPHOSGeometry::ImpactOnEmc(const Double_t theta, const Double_t phi, Int_t & ModuleNumber, Double_t & z, Double_t & x) const
+void AliPHOSGeometry::ImpactOnEmc(Double_t theta, Double_t phi, Int_t & moduleNumber, Double_t & z, Double_t & x) const
 {
   // calculates the impact coordinates on PHOS of a neutral particle  
   // emitted in the direction theta and phi in the ALICE global coordinate system
 
   // searches for the PHOS EMC module
-  ModuleNumber = 0 ; 
+
+  moduleNumber = 0 ; 
   Double_t tm, tM, pm, pM ; 
   Int_t index = 1 ; 
-  while ( ModuleNumber == 0 && index <= GetNModules() ) { 
+  while ( moduleNumber == 0 && index <= GetNModules() ) { 
     EmcModuleCoverage(index, tm, tM, pm, pM) ; 
     if ( (theta >= tm && theta <= tM) && (phi >= pm && phi <= pM ) ) 
-      ModuleNumber = index ; 
+      moduleNumber = index ; 
     index++ ;    
   }
-  if ( ModuleNumber != 0 ) {
-    Float_t phi0 =  GetPHOSAngle(ModuleNumber) *  (TMath::Pi() / 180.) + 1.5 * TMath::Pi()  ;  
+  if ( moduleNumber != 0 ) {
+    Float_t phi0 =  GetPHOSAngle(moduleNumber) *  (TMath::Pi() / 180.) + 1.5 * TMath::Pi()  ;  
     Float_t y0  =  GetIPtoCrystalSurface()  ;   
     Double_t angle = phi - phi0; 
     x = y0 * TMath::Tan(angle) ; 
@@ -343,6 +350,47 @@ void AliPHOSGeometry::ImpactOnEmc(const Double_t theta, const Double_t phi, Int_
   }
 }
 
+//____________________________________________________________________________
+void AliPHOSGeometry::ImpactOnEmc(TVector3 vec, Int_t & moduleNumber, Double_t & z, Double_t & x) const
+{
+  // calculates the impact coordinates on PHOS of a neutral particle  
+  // emitted in the direction theta and phi in the ALICE global coordinate system
+  // searches for the PHOS EMC module
+
+  TParticle p ; 
+  p.SetMomentum(vec.X(), vec.Y(), vec.Z(), 0.) ; 
+  
+  ImpactOnEmc(p, moduleNumber, z, x) ;
+}
+
+//____________________________________________________________________________
+void AliPHOSGeometry::ImpactOnEmc(TParticle p, Int_t & moduleNumber, Double_t & z, Double_t & x) const
+{
+  // calculates the impact coordinates on PHOS of a neutral particle  
+  // emitted in the direction theta and phi in the ALICE global coordinate system
+
+  // searches for the PHOS EMC module
+  Double_t theta = p.Theta() ; 
+  Double_t phi   = p.Phi() ; 
+
+  ImpactOnEmc(theta, phi, moduleNumber, z, x) ;
+}
+
+//____________________________________________________________________________
+Bool_t  AliPHOSGeometry::Impact(const TParticle * particle) const 
+{
+  // Tells if a particle enters PHOS
+  Bool_t in=kFALSE;
+  Int_t moduleNumber=0;
+  Double_t z,x;
+  ImpactOnEmc(particle->Theta(),particle->Phi(),moduleNumber,z,x);
+  if(moduleNumber) 
+    in=kTRUE;
+  else 
+    in=kFALSE;
+  return in;
+}
+
 //____________________________________________________________________________
 Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  AbsId) const
 {
@@ -373,7 +421,7 @@ Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  AbsId) c
 
 //____________________________________________________________________________
 
-void AliPHOSGeometry::RelPosInAlice(const Int_t id, TVector3 & pos ) const
+void AliPHOSGeometry::RelPosInAlice(Int_t id, TVector3 & pos ) const
 {
   // Converts the absolute numbering into the global ALICE coordinate system
   
@@ -399,7 +447,7 @@ void AliPHOSGeometry::RelPosInAlice(const Int_t id, TVector3 & pos ) const
     pos.SetY(y0) ;
     
     Float_t phi           = GetPHOSAngle( phosmodule) ; 
-    Double_t const kRADDEG = 180.0 / kPI ;
+    Double_t const kRADDEG = 180.0 / TMath::Pi() ;
     Float_t rphi          = phi / kRADDEG ; 
     
     TRotation rot ;
@@ -410,6 +458,19 @@ void AliPHOSGeometry::RelPosInAlice(const Int_t id, TVector3 & pos ) const
     pos.Transform(rot) ; // rotate the baby 
 } 
 
+//____________________________________________________________________________
+void AliPHOSGeometry::RelPosToAbsId(const Int_t module, const Double_t x, const Double_t z, Int_t & AbsId) const
+{
+  // converts local PHOS-module (x, z) coordinates to absId 
+  Int_t relid[4] ;
+  relid[0] = module ;
+  relid[1] = 0 ;
+  relid[2] = static_cast<Int_t>(TMath::Ceil( x/ GetCellStep() + GetNPhi() / 2.) );
+  relid[3] = static_cast<Int_t>(TMath::Ceil(-z/ GetCellStep() + GetNZ()   / 2.) ) ;
+  
+  RelToAbsNumbering(relid,AbsId) ;
+}
+
 //____________________________________________________________________________
 void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) const 
 {
@@ -420,12 +481,40 @@ void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t &
   Int_t column     = relid[3] ; //offset along z axis
 
   
-  if ( relid[1] == 0 ) { // its a PbW04 crystal 
-    x = - ( GetNPhi()/2. - row    + 0.5 ) *  GetCrystalSize(0) ; // position of Xtal with respect
-    z =   ( GetNZ()  /2. - column + 0.5 ) *  GetCrystalSize(2) ; // of center of PHOS module  
+  if ( relid[1] == 0 ) { // its a PbW04 crystal
+    x = - ( GetNPhi()/2. - row    + 0.5 ) *  GetCellStep() ; // position of Xtal with respect
+    z = - ( GetNZ()  /2. - column + 0.5 ) *  GetCellStep() ; // of center of PHOS module  
   }  
   else  {    
     x = - ( GetNumberOfCPVPadsPhi()/2. - row    - 0.5 ) * GetPadSizePhi()  ; // position of pad  with respect
-    z =   ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module  
+    z = - ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module  
   }
 }
+
+//____________________________________________________________________________
+
+TVector3 AliPHOSGeometry::GetModuleCenter(char *det, Int_t module) const
+{
+  // Returns a position of the center of the CPV or EMC module
+  Float_t rDet = 0.;
+  if      (strcmp(det,"CPV") == 0) rDet  = GetIPtoCPVDistance   ();
+  else if (strcmp(det,"EMC") == 0) rDet  = GetIPtoCrystalSurface();
+  else Fatal("GetModuleCenter","Wrong detector name %s",det);
+
+  Float_t angle = GetPHOSAngle(module); // (40,20,0,-20,-40) degrees
+  angle *= TMath::Pi()/180;
+  angle += 3*TMath::Pi()/2.;
+  return TVector3(rDet*TMath::Cos(angle), rDet*TMath::Sin(angle), 0.);
+}
+
+//____________________________________________________________________________
+
+TVector3 AliPHOSGeometry::Global2Local(TVector3 globalPosition, Int_t module) const
+{
+  // Transforms a global position of the rec.point to the local coordinate system
+  Float_t angle = GetPHOSAngle(module); // (40,20,0,-20,-40) degrees
+  angle *= TMath::Pi()/180;
+  angle += 3*TMath::Pi()/2.;
+  globalPosition.RotateZ(-angle);
+  return TVector3(globalPosition.Y(),globalPosition.X(),globalPosition.Z());
+}