EMCAL/AliEMCALGeometry.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
  16 /* $Id$*/
  17
  18 //_________________________________________________________________________
  19 // Geometry class  for EMCAL : singleton
  20 // EMCAL consists of layers of scintillator and lead
  21 // Places the the Barrel Geometry of The EMCAL at Midrapidity
  22 // between 0 and 120 degrees of Phi and
  23 // -0.7 to 0.7 in eta
  24 // Number of Modules and Layers may be controlled by
  25 // the name of the instance defined
  26 // EMCALArch2x has more modules along both phi and eta
  27 // EMCALArchxa has less Layers in the Radial Direction
  28 //*-- Author: Sahal Yacoob (LBL / UCT)
  29 //     and  : Yves Schutz (SUBATECH)
  30 //     and  : Jennifer Klay (LBL)
  31
  32 // --- ROOT system ---
  33
  34 // --- Standard library ---
  35
  36 #include <iostream.h>
  37
  38 // --- AliRoot header files ---
  39 #include <TMath.h>
  40 // -- ALICE Headers.
  41 #include "AliConst.h"
  42 // --- EMCAL headers
  43 #include "AliEMCALGeometry.h"
  44
  45 ClassImp(AliEMCALGeometry);
  46
  47 AliEMCALGeometry *AliEMCALGeometry::fgGeom = 0;
  48 Bool_t            AliEMCALGeometry::fgInit = kFALSE;
  49
  50 //______________________________________________________________________
  51 AliEMCALGeometry::~AliEMCALGeometry(void){
  52     // dtor
  53 }
  54 //______________________________________________________________________
  55 void AliEMCALGeometry::Init(void){
  56     // Initializes the EMCAL parameters
  57
  58     fgInit = kFALSE; // Assume failer untill proven otherwise.
  59     if(!(  (strcmp( fName, "EMCALArch1a" ) == 0) |
  60            (strcmp( fName, "EMCALArch1b" ) == 0) |
  61            (strcmp( fName, "EMCALArch2a" ) == 0) |
  62            (strcmp( fName, "EMCALArch2b" ) == 0) )){
  63         cout <<"Instance " << fName << " undefined" << endl;
  64     } // end if
  65     //
  66     if (((strcmp( fName, "EMCALArch1a" ))    == 0) |
  67         ((strcmp( fName, "EMCALArch1b" ))    == 0)){
  68         fNZ         = 96;
  69         fNPhi       = 144;
  70     } // end if
  71     if (((strcmp( fName, "EMCALArch2a" ))    == 0) |
  72         ((strcmp( fName, "EMCALArch2b" ))    == 0)){
  73         fNZ         = 112;
  74         fNPhi       = 168;
  75     } // end if
  76     if (((strcmp( fName, "EMCALArch1a" ))    == 0) |
  77         ((strcmp( fName, "EMCALArch2a" ))    == 0)){
  78         fNLayers    = 21;
  79     } // end if
  80     if (((strcmp( fName, "EMCALArch1b" ))    == 0) |
  81         ((strcmp( fName, "EMCALArch2b" ))    == 0)){
  82         fNLayers    = 25;
  83     } // end if
  84
  85     // geometry
  86     fAirGap         = 5.0; // cm, air gap between EMCAL mother volume and
  87                            // active material.
  88     fAlFrontThick   = 3.18; // cm, Thickness of front Al layer
  89     fPbRadThickness = 0.5; // cm, Thickness of theh Pb radiators.
  90     fPreShowerSintThick = 0.6; // cm, Thickness of the sintilator for the
  91                                // preshower part of the calorimeter
  92     fFullShowerSintThick = 0.5; // cm, Thickness of the sintilator for the
  93                                 // full shower part of the calorimeter
  94     fArm1PhiMin     =  60.0; // degrees, Starting EMCAL Phi position
  95     fArm1PhiMax     = 180.0; // degrees, Ending EMCAL Phi position
  96     fArm1EtaMin     = -0.7; // pseudorapidity, Starting EMCAL Eta position
  97     fArm1EtaMax     = +0.7; // pseudorapidity, Ending EMCAL Eta position
  98     fIPDistance     = 454.0; // cm, Radial distance to inner surface of EMCAL
  99     fShellThickness = GetAlFrontThickness() + 2.*GetPreSintThick() +
 100         (fNLayers-2)*GetFullSintThick()+(fNLayers-1)*GetPbRadThick();
 101     //below; cm, Z lenght of the EMCAL.
 102     fZLength        = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax);
 103     fEnvelop[0]     = fIPDistance; // mother volume inner radius
 104     fEnvelop[1]     = fIPDistance + fShellThickness; // mother volume outer r.
 105     fEnvelop[2]     = 1.00001*fZLength; // add some padding for mother volume.
 106     fGap2Active     = 1.0;  // cm, Gap between
 107     fgInit = kTRUE;
 108 }
 109 //______________________________________________________________________
 110 AliEMCALGeometry *  AliEMCALGeometry::GetInstance(){
 111     // Returns the pointer of the unique instance
 112
 113     return (AliEMCALGeometry *) fgGeom;
 114 }
 115 //______________________________________________________________________
 116 AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
 117                                                 const Text_t* title){
 118     // Returns the pointer of the unique instance
 119
 120     AliEMCALGeometry * rv = 0;
 121     if ( fgGeom == 0 ) {
 122         if ( strcmp(name,"") == 0 ) rv = 0;
 123         else {
 124             fgGeom = new AliEMCALGeometry(name, title);
 125             if ( fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
 126             else {
 127                 rv = 0;
 128                 delete fgGeom;
 129                 fgGeom = 0;
 130             } // end if fgInit
 131         } // end if strcmp(name,"")
 132     }else{
 133         if ( strcmp(fgGeom->GetName(), name) != 0 ) {
 134             cout << "AliEMCALGeometry <E> : current geometry is "
 135                  << fgGeom->GetName() << endl
 136                  << "                      you cannot call     " << name
 137                  << endl;
 138         }else{
 139             rv = (AliEMCALGeometry *) fgGeom;
 140         } // end if
 141     }  // end if fgGeom
 142     return rv;
 143 }
 144 //______________________________________________________________________
 145 Int_t AliEMCALGeometry::TowerIndex(Int_t ieta,Int_t iphi,Int_t ipre) const {
 146     // Returns the tower index number from the based on the Z and Phi
 147     // index numbers. There are 2 times the number of towers to separate
 148     // out the full towsers from the pre-towsers.
 149     // Inputs:
 150     //   Int_t ieta    // index allong z axis [1-fNZ]
 151     //   Int_t iphi  // index allong phi axis [1-fNPhi]
 152     //   Int_t ipre  // 0 = Full tower, 1 = Pre-shower tower only. [0,1]
 153     // Outputs:
 154     //   none.
 155     // Returned
 156     // Int_t the absoulute tower index. [1-2*fNZ*fNPhi]
 157     Int_t index;
 158
 159     if((ieta<=0 || ieta>GetNEta()) || (iphi<=0 || iphi>GetNPhi()) ||
 160        (ipre<0 || ipre>1) ){
 161         cout << "inputs out of range ieta=" << ieta << " [1-" << GetNEta();
 162         cout << "] iphi=" << iphi << " [1-" << GetNPhi() << "] ipre=";
 163         cout << ipre << "[0,1]. returning -1" << endl;
 164         return -1;
 165     } // end if
 166     index = iphi + GetNPhi()*(ieta-1) + ipre*(GetNPhi()*GetNEta());
 167     return index;
 168 }
 169 //______________________________________________________________________
 170 void AliEMCALGeometry::TowerIndexes(Int_t index,Int_t &ieta,Int_t &iphi,
 171                                     Int_t &ipre) const {
 172     // given the tower index number it returns the based on the Z and Phi
 173     // index numbers and if it is for the full tower or the pre-tower number.
 174     // There are 2 times the number of towers to separate
 175     // out the full towsers from the pre-towsers.
 176     // Inputs:
 177     //   Int_t index // Tower index number [1-2*fNZ*fNPhi]
 178     // Outputs:
 179     //   Int_t ieta    // index allong z axis [1-fNZ]
 180     //   Int_t iphi  // index allong phi axis [1-fNPhi]
 181     //   Int_t ipre  // 0 = Full tower, 1 = Pre-shower tower only. [0,1]
 182     // Returned
 183     //   none.
 184     Int_t itowers;
 185
 186     itowers = GetNEta()*GetNPhi();
 187     if(index<1 || index>2*itowers){
 188         cout << "index=" << index <<" is out of range [1-";
 189         cout << 2*itowers << "], returning -1 for all." << endl;
 190         ieta = -1; iphi = -1; ipre = -1;
 191         return ;
 192     } // end if
 193     ipre = 0;
 194     if(index>itowers){ // pre shower indexs
 195         ipre = 1;
 196         index = index - itowers;
 197     } // end if
 198     ieta = 1+ (Int_t)((index-1)/GetNPhi());
 199     iphi = index - GetNPhi()*(ieta-1);
 200     return;
 201 }
 202 //______________________________________________________________________
 203 void AliEMCALGeometry::EtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi) const {
 204     // given the tower index number it returns the based on the eta and phi
 205     // of the tower.
 206     // Inputs:
 207     //   Int_t index // Tower index number [1-2*fNZ*fNPhi]
 208     // Outputs:
 209     //   Float_t eta  // eta of center of tower in pseudorapidity
 210     //   Float_t phi  // phi of center of tower in degrees
 211     // Returned
 212     //   none.
 213     Int_t ieta,iphi,ipre;
 214     Double_t deta,dphi,phid;
 215
 216     TowerIndexes(index,ieta,iphi,ipre);
 217     deta = (GetArm1EtaMax()-GetArm1EtaMin())/((Float_t)GetNEta());
 218     eta  = GetArm1EtaMin() + (((Float_t)ieta)-0.5)*deta;
 219     dphi = (GetArm1PhiMax() - GetArm1PhiMin())/((Float_t)GetNPhi());  // in degrees.
 220     phid = GetArm1PhiMin() + dphi*((Float_t)iphi -0.5);//iphi range [1-fNphi].
 221     phi  = phid;
 222 }
 223 //______________________________________________________________________
 224 Int_t AliEMCALGeometry::TowerIndexFromEtaPhi(Float_t eta,Float_t phi) const {
 225     // returns the tower index number based on the eta and phi of the tower.
 226     // Inputs:
 227     //   Float_t eta  // eta of center of tower in pseudorapidity
 228     //   Float_t phi  // phi of center of tower in degrees
 229     // Outputs:
 230     //   none.
 231     // Returned
 232     //   Int_t index // Tower index number [1-fNZ*fNPhi]
 233     Int_t ieta,iphi;
 234
 235     ieta = 1 + (Int_t)(((Float_t)GetNEta())*(eta-GetArm1EtaMin())/
 236                   (GetArm1EtaMax() - GetArm1EtaMin()));
 237     if(ieta<=0 || ieta>GetNEta()){
 238         cout << "TowerIndexFromEtaPhi:";
 239         cout << "ieta = "<< ieta << " eta=" << eta << " is outside of EMCAL. etamin=";
 240         cout << GetArm1EtaMin() << " to etamax=" << GetArm1EtaMax();
 241         cout << " returning -1" << endl;
 242         return -1;
 243     } // end if
 244     iphi = 1 + (Int_t)(((Float_t)GetNPhi())*(phi-GetArm1PhiMin())/
 245                   ((Float_t)(GetArm1PhiMax() - GetArm1PhiMin())));
 246     if(iphi<=0 || iphi>GetNPhi()){
 247         cout << "TowerIndexFromEtaPhi:";
 248         cout << "iphi=" << iphi << " phi=" << phi << " is outside of EMCAL.";
 249         cout << " Phimin=" << GetArm1PhiMin() << " PhiMax=" << GetArm1PhiMax();
 250         cout << " returning -1" << endl;
 251         return -1;
 252     } // end if
 253     return TowerIndex(ieta,iphi,0);
 254 }
 255 //______________________________________________________________________
 256 Int_t AliEMCALGeometry::PreTowerIndexFromEtaPhi(Float_t eta,Float_t phi) const {
 257     // returns the pretower index number based on the eta and phi of the tower.
 258     // Inputs:
 259     //   Float_t eta  // eta of center of tower in pseudorapidity
 260     //   Float_t phi  // phi of center of tower in degrees
 261     // Outputs:
 262     //   none.
 263     // Returned
 264     //   Int_t index // PreTower index number [fNZ*fNPhi-2*fNZ*fNPhi]
 265
 266     return GetNEta()*GetNPhi()+TowerIndexFromEtaPhi(eta,phi);
 267 }
 268 //______________________________________________________________________
 269 Bool_t AliEMCALGeometry::AbsToRelNumbering(Int_t AbsId, Int_t *relid) const {
 270     // Converts the absolute numbering into the following array/
 271     //  relid[0] = EMCAL Arm number 1:1
 272     //  relid[1] = 0  Not in Pre Shower layers
 273     //           = -1 In Pre Shower
 274     //  relid[2] = Row number inside EMCAL
 275     //  relid[3] = Column number inside EMCAL
 276     // Input:
 277     //   Int_t AbsId // Tower index number [1-2*fNZ*fNPhi]
 278     // Outputs:
 279     //   Int_t *relid // array of 5. Discribed above.
 280     Bool_t rv  = kTRUE ;
 281     Int_t ieta=0,iphi=0,ipre=0,index=AbsId;
 282
 283     TowerIndexes(index,ieta,iphi,ipre);
 284     relid[0] = 1;
 285     relid[1] = 0;
 286     if(ipre==1)
 287       relid[1] = -1;
 288     relid[2] = ieta;
 289     relid[3] = iphi;
 290
 291     return rv;
 292 }
 293 //______________________________________________________________________
 294 void AliEMCALGeometry::PosInAlice(const Int_t *relid,Float_t &theta,
 295                                      Float_t &phi) const {
 296     // Converts the relative numbering into the local EMCAL-module (x, z)
 297     // coordinates
 298     Int_t ieta   = relid[2]; // offset along x axis
 299     Int_t iphi = relid[3]; // offset along z axis
 300     Int_t ipre = relid[1]; // indicates -1 preshower, or 0 full tower.
 301     Int_t index;
 302     Float_t eta;
 303
 304     if(ipre==-1) ipre = 1;
 305     index = TowerIndex(ieta,iphi,ipre);
 306     EtaPhiFromIndex(index,eta,phi);
 307     theta = 180.*(2.0*TMath::ATan(TMath::Exp(-eta)))/TMath::Pi();
 308
 309     return;
 310 }
 311
 312 //______________________________________________________________________
 313 void AliEMCALGeometry::XYZFromIndex(const Int_t *relid,Float_t &x,Float_t &y, Float_t &z) const {
 314     // given the tower relative number it returns the X, Y and Z
 315     // of the tower.
 316
 317     // Outputs:
 318     //   Float_t x  // x of center of tower in cm
 319     //   Float_t y  // y of center of tower in cm
 320     //   Float_t z  // z of centre of tower in cm
 321     // Returned
 322     //   none.
 323
 324     Float_t eta,theta, phi,cyl_radius,kDeg2Rad;
 325
 326     Int_t ieta   = relid[2]; // offset along x axis
 327     Int_t iphi = relid[3]; // offset along z axis
 328     Int_t ipre = relid[1]; // indicates -1 preshower, or 0 full tower.
 329     Int_t index;
 330
 331
 332     if(ipre==-1) ipre = 1;
 333     index = TowerIndex(ieta,iphi,ipre);
 334     EtaPhiFromIndex(index,eta,phi);
 335     theta = 180.*(2.0*TMath::ATan(TMath::Exp(-eta)))/TMath::Pi();
 336
 337
 338
 339     kDeg2Rad = TMath::Pi() / static_cast<Double_t>(180) ;
 340     cyl_radius = GetIPDistance()+ GetAirGap() ;
 341     x =  cyl_radius * TMath::Cos(phi * kDeg2Rad ) ;
 342     y =  cyl_radius * TMath::Cos(phi * kDeg2Rad ) ;
 343     z =  cyl_radius / TMath::Tan(theta * kDeg2Rad ) ;
 344
 345  return;
 346 }
 347
 348 //______________________________________________________________________
 349 /*
 350 Boot_t AliEMCALGeometry::AreNeighbours(Int_t index1,Int_t index2) const {
 351     // Returns kTRUE if the two towers are neighbours or not, including
 352     // diagonals. Both indexes are required to be either towers or preshower.
 353     // Inputs:
 354     //   Int_t index1  // index of tower 1
 355     //   Int_t index2  // index of tower 2
 356     // Outputs:
 357     //   none.
 358     // Returned
 359     //   Boot_t kTRUE if the towers are neighbours otherwise false.
 360     Boot_t anb = kFALSE;
 361     Int_t ieta1 = 0, ieta2 = 0, iphi1 = 0, iphi2 = 0, ipre1 = 0, ipre2 = 0;
 362
 363     TowerIndexes(index1,ieta1,iphi1,ipre1);
 364     TowerIndexes(index2,ieta2,iphi2,ipre2);
 365     if(ipre1!=ipre2) return anb;
 366     if((ieta1>=ieta2-1 && ieta1<=ieta2+1) && (iphi1>=iphi2-1 &&iphi1<=iphi2+1))
 367                                                                  anb = kTRUE;
 368     return anb;
 369 }
 370  */