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 80 and 180(or 190) 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 //     EMCAL geometry tree:
  27 //     EMCAL -> superModule -> module -> cell
  28 //     Indexes
  29 //     absId -> nSupMod     -> nTower -> (nIphi,nIeta)
  30 //
  31 //*-- Author: Sahal Yacoob (LBL / UCT)
  32 //     and  : Yves Schutz (SUBATECH)
  33 //     and  : Jennifer Klay (LBL)
  34 //     SHASHLYK : Aleksei Pavlinov (WSU)
  35 //
  36 // --- AliRoot header files ---
  37 #include <assert.h>
  38 #include <Riostream.h>
  39
  40 #include <TMath.h>
  41 #include <TVector3.h>
  42 #include <TObjArray.h>
  43 #include <TGeoManager.h>
  44 #include <TGeoNode.h>
  45 #include <TGeoMatrix.h>
  46 #include <TMatrixD.h>
  47 #include <TObjString.h>
  48 #include <TClonesArray.h>
  49 #include <TBrowser.h>
  50
  51 // -- ALICE Headers.
  52 //#include "AliConst.h"
  53 #include "AliLog.h"
  54
  55 // --- EMCAL headers
  56 #include "AliEMCALGeometry.h"
  57 #include "AliEMCALShishKebabTrd1Module.h"
  58 #include "AliEMCALRecPoint.h"
  59 #include "AliEMCALDigit.h"
  60 #include "AliEMCALHistoUtilities.h"
  61
  62 ClassImp(AliEMCALGeometry)
  63
  64 // these initialisations are needed for a singleton
  65 AliEMCALGeometry  *AliEMCALGeometry::fgGeom      = 0;
  66 Bool_t             AliEMCALGeometry::fgInit      = kFALSE;
  67
  68
  69 AliEMCALGeometry::AliEMCALGeometry()
  70   : AliGeometry(),
  71     fGeoName(0),fArrayOpts(0),fAlFrontThick(0.),fECPbRadThickness(0.),fECScintThick(0.),
  72     fNECLayers(0),fArm1PhiMin(0.),fArm1PhiMax(0.),fArm1EtaMin(0.),fArm1EtaMax(0.),fIPDistance(0.),
  73     fShellThickness(0.),fZLength(0.),fGap2Active(0.),fNZ(0),fNPhi(0),fSampling(0.),fNumberOfSuperModules(0),
  74     fSteelFrontThick(0.),fFrontSteelStrip(0.),fLateralSteelStrip(0.),fPassiveScintThick(0.),fPhiModuleSize(0.),
  75     fEtaModuleSize(0.),fPhiTileSize(0.),fEtaTileSize(0.),fLongModuleSize(0.),fNPhiSuperModule(0),fNPHIdiv(0),fNETAdiv(0),
  76     fNCells(0),fNCellsInSupMod(0),fNCellsInTower(0),fNTRU(0),fNTRUEta(0),fNTRUPhi(0),fTrd1Angle(0.),f2Trd1Dx2(0.),
  77     fPhiGapForSM(0.),fKey110DEG(0),fPhiBoundariesOfSM(0), fPhiCentersOfSM(0),fEtaMaxOfTRD1(0),
  78     fTrd2AngleY(0.),f2Trd2Dy2(0.),fEmptySpace(0.),fTubsR(0.),fTubsTurnAngle(0.),fCentersOfCellsEtaDir(0),
  79     fCentersOfCellsXDir(0),fCentersOfCellsPhiDir(0),fEtaCentersOfCells(0),fPhiCentersOfCells(0),
  80     fShishKebabTrd1Modules(0),fNAdditionalOpts(0)
  81 {
  82   // default ctor only for internal usage (singleton)
  83   // must be kept public for root persistency purposes, but should never be called by the outside world
  84   //  CreateListOfTrd1Modules();
  85   AliDebug(2, "AliEMCALGeometry : default ctor ");
  86 }
  87 //______________________________________________________________________
  88 AliEMCALGeometry::AliEMCALGeometry(const Text_t* name, const Text_t* title)
  89   : AliGeometry(name, title),
  90     fGeoName(0),fArrayOpts(0),fAlFrontThick(0.),fECPbRadThickness(0.),fECScintThick(0.),
  91     fNECLayers(0),fArm1PhiMin(0.),fArm1PhiMax(0.),fArm1EtaMin(0.),fArm1EtaMax(0.),fIPDistance(0.),
  92     fShellThickness(0.),fZLength(0.),fGap2Active(0.),fNZ(0),fNPhi(0),fSampling(0.),fNumberOfSuperModules(0),
  93     fSteelFrontThick(0.),fFrontSteelStrip(0.),fLateralSteelStrip(0.),fPassiveScintThick(0.),fPhiModuleSize(0.),
  94     fEtaModuleSize(0.),fPhiTileSize(0.),fEtaTileSize(0.),fLongModuleSize(0.),fNPhiSuperModule(0),fNPHIdiv(0),fNETAdiv(0),
  95     fNCells(0),fNCellsInSupMod(0),fNCellsInTower(0),fNTRU(0),fNTRUEta(0),fNTRUPhi(0),fTrd1Angle(0.),f2Trd1Dx2(0.),
  96     fPhiGapForSM(0.),fKey110DEG(0),fPhiBoundariesOfSM(0), fPhiCentersOfSM(0), fEtaMaxOfTRD1(0),
  97     fTrd2AngleY(0.),f2Trd2Dy2(0.),fEmptySpace(0.),fTubsR(0.),fTubsTurnAngle(0.),fCentersOfCellsEtaDir(0),
  98     fCentersOfCellsXDir(0),fCentersOfCellsPhiDir(0),fEtaCentersOfCells(0),fPhiCentersOfCells(0),
  99     fShishKebabTrd1Modules(0),fNAdditionalOpts(0)
 100 {
 101   // ctor only for internal usage (singleton)
 102   AliDebug(2, Form("AliEMCALGeometry(%s,%s) ", name,title));
 103
 104   Init();
 105
 106   CreateListOfTrd1Modules();
 107
 108   if (AliDebugLevel()>=2) {
 109     PrintGeometry();
 110   }
 111
 112 }
 113 //______________________________________________________________________
 114 AliEMCALGeometry::AliEMCALGeometry(const AliEMCALGeometry& geom)
 115   : AliGeometry(geom),
 116     fGeoName(geom.fGeoName),
 117     fArrayOpts(geom.fArrayOpts),
 118     fAlFrontThick(geom.fAlFrontThick),
 119     fECPbRadThickness(geom.fECPbRadThickness),
 120     fECScintThick(geom.fECScintThick),
 121     fNECLayers(geom.fNECLayers),
 122     fArm1PhiMin(geom.fArm1PhiMin),
 123     fArm1PhiMax(geom.fArm1PhiMax),
 124     fArm1EtaMin(geom.fArm1EtaMin),
 125     fArm1EtaMax(geom.fArm1EtaMax),
 126     fIPDistance(geom.fIPDistance),
 127     fShellThickness(geom.fShellThickness),
 128     fZLength(geom.fZLength),
 129     fGap2Active(geom.fGap2Active),
 130     fNZ(geom.fNZ),
 131     fNPhi(geom.fNPhi),
 132     fSampling(geom.fSampling),
 133     fNumberOfSuperModules(geom.fNumberOfSuperModules),
 134     fSteelFrontThick(geom.fSteelFrontThick),
 135     fFrontSteelStrip(geom.fFrontSteelStrip),
 136     fLateralSteelStrip(geom.fLateralSteelStrip),
 137     fPassiveScintThick(geom.fPassiveScintThick),
 138     fPhiModuleSize(geom.fPhiModuleSize),
 139     fEtaModuleSize(geom.fEtaModuleSize),
 140     fPhiTileSize(geom.fPhiTileSize),
 141     fEtaTileSize(geom.fEtaTileSize),
 142     fLongModuleSize(geom.fLongModuleSize),
 143     fNPhiSuperModule(geom.fNPhiSuperModule),
 144     fNPHIdiv(geom.fNPHIdiv),
 145     fNETAdiv(geom.fNETAdiv),
 146     fNCells(geom.fNCells),
 147     fNCellsInSupMod(geom.fNCellsInSupMod),
 148     fNCellsInTower(geom.fNCellsInTower),
 149     fNTRU(geom.fNTRU),
 150     fNTRUEta(geom.fNTRUEta),
 151     fNTRUPhi(geom.fNTRUPhi),
 152     fTrd1Angle(geom.fTrd1Angle),
 153     f2Trd1Dx2(geom.f2Trd1Dx2),
 154     fPhiGapForSM(geom.fPhiGapForSM),
 155     fKey110DEG(geom.fKey110DEG),
 156     fPhiBoundariesOfSM(geom.fPhiBoundariesOfSM),
 157     fPhiCentersOfSM(geom.fPhiCentersOfSM),
 158     fEtaMaxOfTRD1(geom.fEtaMaxOfTRD1),
 159     fTrd2AngleY(geom.fTrd2AngleY),
 160     f2Trd2Dy2(geom.f2Trd2Dy2),
 161     fEmptySpace(geom.fEmptySpace),
 162     fTubsR(geom.fTubsR),
 163     fTubsTurnAngle(geom.fTubsTurnAngle),
 164     fCentersOfCellsEtaDir(geom.fCentersOfCellsEtaDir),
 165     fCentersOfCellsXDir(geom.fCentersOfCellsXDir),
 166     fCentersOfCellsPhiDir(geom.fCentersOfCellsPhiDir),
 167     fEtaCentersOfCells(geom.fEtaCentersOfCells),
 168     fPhiCentersOfCells(geom.fPhiCentersOfCells),
 169     fShishKebabTrd1Modules(geom.fShishKebabTrd1Modules),
 170     fNAdditionalOpts(geom.fNAdditionalOpts)
 171 {
 172   //copy ctor
 173 }
 174
 175 //______________________________________________________________________
 176 AliEMCALGeometry::~AliEMCALGeometry(void){
 177     // dtor
 178 }
 179 //______________________________________________________________________
 180 void AliEMCALGeometry::Init(void){
 181   // Initializes the EMCAL parameters
 182   // naming convention : GUV_WX_N_ gives the composition of a tower
 183   // WX inform about the composition of the EM calorimeter section:
 184   //   thickness in mm of Pb radiator (W) and of scintillator (X), and number of scintillator layers (N)
 185   // New geometry: EMCAL_55_25
 186   // 24-aug-04 for shish-kebab
 187   // SHISH_25 or SHISH_62
 188   // 11-oct-05   - correction for pre final design
 189   // Feb 06,2006 - decrease the weight of EMCAL
 190   //
 191   // Oct 30,2006 - SHISH_TRD1_CURRENT_1X1, SHISH_TRD1_CURRENT_2X2 or SHISH_TRD1_CURRENT_3X3;
 192   //
 193
 194   fAdditionalOpts[0] = "nl=";    // number of sampling layers (fNECLayers)
 195   fAdditionalOpts[1] = "pbTh=";  // cm, Thickness of the Pb   (fECPbRadThick)
 196   fAdditionalOpts[2] = "scTh=";  // cm, Thickness of the Sc    (fECScintThick)
 197   fAdditionalOpts[3] = "latSS=";  // cm, Thickness of lateral steel strip (fLateralSteelStrip)
 198
 199   fNAdditionalOpts = sizeof(fAdditionalOpts) / sizeof(char*);
 200
 201   fgInit = kFALSE; // Assume failed until proven otherwise.
 202   fGeoName   = GetName();
 203   fGeoName.ToUpper();
 204   fKey110DEG = 0;
 205   if(fGeoName.Contains("110DEG") || fGeoName.Contains("CURRENT")) fKey110DEG = 1; // for GetAbsCellId
 206   fShishKebabTrd1Modules = 0;
 207   fTrd2AngleY = f2Trd2Dy2 = fEmptySpace = fTubsR = fTubsTurnAngle = 0;
 208
 209   fNZ             = 114;        // granularity along Z (eta)
 210   fNPhi           = 168;        // granularity in phi (azimuth)
 211   fArm1PhiMin     = 80.0;       // degrees, Starting EMCAL Phi position
 212   fArm1PhiMax     = 190.0;      // degrees, Ending EMCAL Phi position
 213   fArm1EtaMin     = -0.7;       // pseudorapidity, Starting EMCAL Eta position
 214   fArm1EtaMax     = +0.7;       // pseudorapidity, Ending EMCAL Eta position
 215   fIPDistance     = 454.0;      // cm, Radial distance to inner surface of EMCAL
 216   fPhiGapForSM    = 0.;         // cm, only for final TRD1 geometry
 217   for(int i=0; i<12; i++) fMatrixOfSM[i] = 0;
 218
 219   // geometry
 220   if(fGeoName.Contains("SHISH")){ // Only shahslyk now
 221     // 7-sep-05; integration issue
 222     fArm1PhiMin     = 80.0;     // 60  -> 80
 223     fArm1PhiMax     = 180.0;    // 180 -> 190
 224
 225     fNumberOfSuperModules = 10; // 12 = 6 * 2 (6 in phi, 2 in Z);
 226     fSteelFrontThick = 2.54;    //  9-sep-04
 227     fIPDistance      = 460.0;
 228     fFrontSteelStrip = fPassiveScintThick = 0.0; // 13-may-05
 229     fLateralSteelStrip = 0.025; // before MAY 2005
 230     fPhiModuleSize   = fEtaModuleSize   = 11.4;
 231     fPhiTileSize = fEtaTileSize      = 5.52; // (11.4-5.52*2)/2. = 0.18 cm (wall thickness)
 232     fNPhi            = 14;
 233     fNZ              = 30;
 234     fAlFrontThick    = fGap2Active = 0;
 235     fNPHIdiv = fNETAdiv = 2;
 236
 237     fNECLayers       = 62;
 238     fECScintThick    = fECPbRadThickness = 0.2;
 239     fSampling        = 1.;  // 30-aug-04 - should be calculated
 240     if(fGeoName.Contains("TWIST")) { // all about EMCAL module
 241       fNZ             = 27;  // 16-sep-04
 242     } else if(fGeoName.Contains("TRD")) {
 243       fIPDistance      = 428.0;  //  11-may-05
 244       fSteelFrontThick = 0.0;    // 3.17 -> 0.0; 28-mar-05 : no stell plate
 245       fNPhi            = 12;
 246       fSampling       = 12.327;
 247       fPhiModuleSize = fEtaModuleSize = 12.26;
 248       fNZ            = 26;     // 11-oct-04
 249       fTrd1Angle     = 1.3;    // in degree
 250 // 18-nov-04; 1./0.08112=12.327
 251 // http://pdsfweb01.nersc.gov/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
 252       if(fGeoName.Contains("TRD1")) {       // 30-jan-05
 253         // for final design
 254         fPhiGapForSM    = 2.;         // cm, only for final TRD1 geometry
 255         if(fGeoName.Contains("MAY05") || fGeoName.Contains("WSUC") || fGeoName.Contains("FINAL") || fGeoName.Contains("CURRENT")){
 256           fNumberOfSuperModules = 12; // 20-may-05
 257           if(fGeoName.Contains("WSUC")) fNumberOfSuperModules = 1; // 27-may-05
 258           fNECLayers     = 77;       // (13-may-05 from V.Petrov)
 259           fPhiModuleSize = 12.5;     // 20-may-05 - rectangular shape
 260           fEtaModuleSize = 11.9;
 261           fECScintThick  = fECPbRadThickness = 0.16;// (13-may-05 from V.Petrov)
 262           fFrontSteelStrip   = 0.025;// 0.025cm = 0.25mm  (13-may-05 from V.Petrov)
 263           fLateralSteelStrip = 0.01; // 0.01cm  = 0.1mm   (13-may-05 from V.Petrov) - was 0.025
 264           fPassiveScintThick = 0.8;  // 0.8cm   = 8mm     (13-may-05 from V.Petrov)
 265           fNZ                = 24;
 266           fTrd1Angle         = 1.5;  // 1.3 or 1.5
 267
 268           if(fGeoName.Contains("FINAL") || fGeoName.Contains("CURRENT")) { // 9-sep-05
 269             fNumberOfSuperModules = 10;
 270             if(GetKey110DEG()) {
 271               fNumberOfSuperModules = 12;// last two modules have size 10 degree in phi (180<phi<190)
 272               fArm1PhiMax = 200.0;       // for XEN1 and turn angle of super modules
 273             }
 274             if(fGeoName.Contains("FINAL")) {
 275               fPhiModuleSize = 12.26 - fPhiGapForSM / Float_t(fNPhi); // first assumption
 276             } else if(fGeoName.Contains("CURRENT")) {
 277               fECScintThick      = 0.176; // 10% of weight reduction
 278               fECPbRadThickness  = 0.144; //
 279               fLateralSteelStrip = 0.015; // 0.015cm  = 0.15mm (Oct 30, from Fred)
 280               fPhiModuleSize     = 12.00;
 281               fPhiGapForSM       = (12.26 - fPhiModuleSize)*fNPhi; // have to check
 282             }
 283             fEtaModuleSize = fPhiModuleSize;
 284             if(fGeoName.Contains("HUGE")) fNECLayers *= 3; // 28-oct-05 for analysing leakage
 285           }
 286         }
 287       } else if(fGeoName.Contains("TRD2")) {       // 30-jan-05
 288         fSteelFrontThick = 0.0;         // 11-mar-05
 289         fIPDistance+= fSteelFrontThick; // 1-feb-05 - compensate absence of steel plate
 290         fTrd1Angle  = 1.64;             // 1.3->1.64
 291         fTrd2AngleY = fTrd1Angle;       //  symmetric case now
 292         fEmptySpace    = 0.2; // 2 mm
 293         fTubsR         = fIPDistance; // 31-jan-05 - as for Fred case
 294
 295         fPhiModuleSize  = fTubsR*2.*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.);
 296         fPhiModuleSize -= fEmptySpace/2.; // 11-mar-05
 297         fEtaModuleSize  = fPhiModuleSize; // 20-may-05
 298         fTubsTurnAngle  = 3.;
 299       }
 300       fNPHIdiv = fNETAdiv  = 2;   // 13-oct-04 - division again
 301       if(fGeoName.Contains("3X3")) {   // 23-nov-04
 302         fNPHIdiv = fNETAdiv  = 3;
 303       } else if(fGeoName.Contains("4X4")) {
 304         fNPHIdiv = fNETAdiv  = 4;
 305       } else if(fGeoName.Contains("1X1")) {
 306         fNPHIdiv = fNETAdiv  = 1;
 307       }
 308     }
 309     if(fGeoName.Contains("25")){
 310       fNECLayers     = 25;
 311       fECScintThick  = fECPbRadThickness = 0.5;
 312     }
 313     if(fGeoName.Contains("WSUC")){ // 18-may-05 - about common structure
 314       fShellThickness = 30.; // should be change
 315       fNPhi = fNZ = 4;
 316     }
 317
 318     CheckAdditionalOptions();
 319     DefineSamplingFraction();
 320
 321     fPhiTileSize = fPhiModuleSize/double(fNPHIdiv) - fLateralSteelStrip; // 13-may-05
 322     fEtaTileSize = fEtaModuleSize/double(fNETAdiv) - fLateralSteelStrip; // 13-may-05
 323
 324     // constant for transition absid <--> indexes
 325     fNCellsInTower  = fNPHIdiv*fNETAdiv;
 326     fNCellsInSupMod = fNCellsInTower*fNPhi*fNZ;
 327     fNCells         = fNCellsInSupMod*fNumberOfSuperModules;
 328     if(GetKey110DEG()) fNCells -= fNCellsInSupMod;
 329
 330     fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness);
 331     if(fGeoName.Contains("MAY05")) fLongModuleSize += (fFrontSteelStrip + fPassiveScintThick);
 332
 333     // 30-sep-04
 334     if(fGeoName.Contains("TRD")) {
 335       f2Trd1Dx2 = fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.);
 336       if(fGeoName.Contains("TRD2")) {  // 27-jan-05
 337         f2Trd2Dy2 = fPhiModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.);
 338       }
 339     }
 340   } else Fatal("Init", "%s is an undefined geometry!", fGeoName.Data()) ;
 341
 342   fNPhiSuperModule = fNumberOfSuperModules/2;
 343   if(fNPhiSuperModule<1) fNPhiSuperModule = 1;
 344
 345   fShellThickness = fAlFrontThick + fGap2Active + fNECLayers*GetECScintThick()+(fNECLayers-1)*GetECPbRadThick();
 346   if(fGeoName.Contains("SHISH")) {
 347     fShellThickness = fSteelFrontThick + fLongModuleSize;
 348     if(fGeoName.Contains("TWIST")) { // 13-sep-04
 349       fShellThickness  = TMath::Sqrt(fLongModuleSize*fLongModuleSize + fPhiModuleSize*fEtaModuleSize);
 350       fShellThickness += fSteelFrontThick;
 351     } else if(fGeoName.Contains("TRD")) { // 1-oct-04
 352       fShellThickness  = TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2);
 353       fShellThickness += fSteelFrontThick;
 354       // Local coordinates
 355       fParSM[0] = GetShellThickness()/2.;
 356       fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
 357       fParSM[2] = 350./2.;
 358     }
 359   }
 360
 361   fZLength        = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax); // Z coverage
 362   fEnvelop[0]     = fIPDistance; // mother volume inner radius
 363   fEnvelop[1]     = fIPDistance + fShellThickness; // mother volume outer r.
 364   fEnvelop[2]     = 1.00001*fZLength; // add some padding for mother volume.
 365
 366   fNumberOfSuperModules = 12;
 367
 368   // SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006
 369   fPhiBoundariesOfSM.Set(fNumberOfSuperModules);
 370   fPhiCentersOfSM.Set(fNumberOfSuperModules/2);
 371   fPhiBoundariesOfSM[0] = TMath::PiOver2() - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules)
 372   fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance);
 373   fPhiCentersOfSM[0]     = TMath::PiOver2();
 374   for(int i=1; i<=4; i++) { // from 2th ro 9th
 375     fPhiBoundariesOfSM[2*i]   = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i;
 376     fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i;
 377     fPhiCentersOfSM[i]         = fPhiCentersOfSM[0]     + 20.*TMath::DegToRad()*i;
 378   }
 379   fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad();
 380   fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance);
 381   fPhiCentersOfSM[5]      = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.;
 382
 383   fgInit = kTRUE;
 384
 385   //TRU parameters. These parameters values are not the final ones.
 386   fNTRU    = 3 ;
 387   fNTRUEta = 3 ;
 388   fNTRUPhi = 1 ;
 389
 390 }
 391
 392 void AliEMCALGeometry::PrintGeometry()
 393 {
 394   // Separate routine is callable from broswer; Nov 7,2006
 395   printf("\nInit: geometry of EMCAL named %s is as follows:\n", fGeoName.Data());
 396   printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
 397   printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f -> for EMCAL envelope only\n",
 398            GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() );
 399
 400   printf( "               ECAL      : %d x (%f cm Pb, %f cm Sc) \n",
 401   GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ;
 402   printf("                fSampling %5.2f \n",  fSampling );
 403   if(fGeoName.Contains("SHISH")){
 404     printf(" fIPDistance       %6.3f cm \n", fIPDistance);
 405     if(fSteelFrontThick>0.)
 406     printf(" fSteelFrontThick  %6.3f cm \n", fSteelFrontThick);
 407     printf(" fNPhi %i   |  fNZ %i \n", fNPhi, fNZ);
 408     printf(" fNCellsInTower %i : fNCellsInSupMod %i : fNCells %i\n",fNCellsInTower, fNCellsInSupMod, fNCells);
 409     if(fGeoName.Contains("MAY05")){
 410       printf(" fFrontSteelStrip         %6.4f cm (thickness of front steel strip)\n",
 411       fFrontSteelStrip);
 412       printf(" fLateralSteelStrip       %6.4f cm (thickness of lateral steel strip)\n",
 413       fLateralSteelStrip);
 414       printf(" fPassiveScintThick  %6.4f cm (thickness of front passive Sc tile)\n",
 415       fPassiveScintThick);
 416     }
 417     printf(" X:Y module size     %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize);
 418     printf(" X:Y   tile size     %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize);
 419     printf(" #of sampling layers %i(fNECLayers) \n", fNECLayers);
 420     printf(" fLongModuleSize     %6.3f cm \n", fLongModuleSize);
 421     printf(" #supermodule in phi direction %i \n", fNPhiSuperModule );
 422   }
 423   if(fGeoName.Contains("TRD")) {
 424     printf(" fTrd1Angle %7.4f\n", fTrd1Angle);
 425     printf(" f2Trd1Dx2  %7.4f\n",  f2Trd1Dx2);
 426     if(fGeoName.Contains("TRD2")) {
 427       printf(" fTrd2AngleY     %7.4f\n", fTrd2AngleY);
 428       printf(" f2Trd2Dy2       %7.4f\n", f2Trd2Dy2);
 429       printf(" fTubsR          %7.2f cm\n", fTubsR);
 430       printf(" fTubsTurnAngle  %7.4f\n", fTubsTurnAngle);
 431       printf(" fEmptySpace     %7.4f cm\n", fEmptySpace);
 432     } else if(fGeoName.Contains("TRD1")){
 433       printf("SM dimensions(TRD1) : dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n",
 434       fParSM[0],fParSM[1],fParSM[2]);
 435       printf(" fPhiGapForSM  %7.4f cm (%7.4f <- phi size in degree)\n",
 436       fPhiGapForSM, TMath::ATan2(fPhiGapForSM,fIPDistance)*TMath::RadToDeg());
 437       if(GetKey110DEG()) printf(" Last two modules have size 10 degree in  phi (180<phi<190)\n");
 438       printf(" phi SM boundaries \n");
 439       for(int i=0; i<fPhiBoundariesOfSM.GetSize()/2.; i++) {
 440         printf(" %i : %7.5f(%7.2f) -> %7.5f(%7.2f) : center %7.5f(%7.2f) \n", i,
 441         fPhiBoundariesOfSM[2*i], fPhiBoundariesOfSM[2*i]*TMath::RadToDeg(),
 442                fPhiBoundariesOfSM[2*i+1], fPhiBoundariesOfSM[2*i+1]*TMath::RadToDeg(),
 443                fPhiCentersOfSM[i], fPhiCentersOfSM[i]*TMath::RadToDeg());
 444       }
 445       printf(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n",
 446                fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1);
 447
 448       printf("\n Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize());
 449       for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
 450         printf(" ind %2.2i : z %8.3f : x %8.3f \n", i,
 451                fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i));
 452         int ind=0; // Nov 21,2006
 453         for(Int_t iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
 454           ind = iphi*fCentersOfCellsEtaDir.GetSize() + i;
 455           printf("%6.4f ", fEtaCentersOfCells[ind]);
 456           if((iphi+1)%12 == 0) printf("\n");
 457         }
 458         printf("\n");
 459       }
 460
 461       printf("\n Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize());
 462       for(Int_t i=0; i<fCentersOfCellsPhiDir.GetSize(); i++) {
 463         double phi=fPhiCentersOfCells.At(i);
 464         printf(" ind %2.2i : y %8.3f : phi %7.5f(%6.2f) \n", i, fCentersOfCellsPhiDir.At(i),
 465                phi, phi*TMath::RadToDeg());
 466       }
 467     }
 468   }
 469 }
 470
 471 void AliEMCALGeometry::PrintCellIndexes(Int_t absId, int pri, char *tit)
 472 {
 473   // Service methods
 474   Int_t nSupMod, nTower, nIphi, nIeta;
 475   Int_t iphi, ieta;
 476   TVector3 vg;
 477
 478   GetCellIndex(absId,  nSupMod, nTower, nIphi, nIeta);
 479   printf(" %s | absId : %i -> nSupMod %i nTower %i nIphi %i nIeta %i \n", tit, absId,  nSupMod, nTower, nIphi, nIeta);
 480   if(pri>0) {
 481     GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
 482     printf(" local SM index : iphi %i : ieta %i \n", iphi,ieta);
 483     GetGlobal(absId, vg);
 484     printf(" vglob : mag %7.2f : perp %7.2f : z %7.2f : eta %6.4f : phi %6.4f(%6.2f) \n",
 485            vg.Mag(), vg.Perp(), vg.Z(), vg.Eta(), vg.Phi(), vg.Phi()*TMath::RadToDeg());
 486   }
 487 }
 488
 489 //______________________________________________________________________
 490 void AliEMCALGeometry::CheckAdditionalOptions()
 491 {
 492   // Feb 06,2006
 493   //Additional options that
 494   //can be used to select
 495   //the specific geometry of
 496   //EMCAL to run
 497
 498   fArrayOpts = new TObjArray;
 499   Int_t nopt = AliEMCALHistoUtilities::ParseString(fGeoName, *fArrayOpts);
 500   if(nopt==1) { // no aditional option(s)
 501     fArrayOpts->Delete();
 502     delete fArrayOpts;
 503     fArrayOpts = 0;
 504     return;
 505   }
 506   for(Int_t i=1; i<nopt; i++){
 507     TObjString *o = (TObjString*)fArrayOpts->At(i);
 508
 509     TString addOpt = o->String();
 510     Int_t indj=-1;
 511     for(Int_t j=0; j<fNAdditionalOpts; j++) {
 512       TString opt = fAdditionalOpts[j];
 513       if(addOpt.Contains(opt,TString::kIgnoreCase)) {
 514           indj = j;
 515         break;
 516       }
 517     }
 518     if(indj<0) {
 519       AliDebug(2,Form("<E> option |%s| unavailable : ** look to the file AliEMCALGeometry.h **\n",
 520                       addOpt.Data()));
 521       assert(0);
 522     } else {
 523       AliDebug(2,Form("<I> option |%s| is valid : number %i : |%s|\n",
 524                       addOpt.Data(), indj, fAdditionalOpts[indj]));
 525       if       (addOpt.Contains("NL=",TString::kIgnoreCase))   {// number of sampling layers
 526         sscanf(addOpt.Data(),"NL=%i", &fNECLayers);
 527         AliDebug(2,Form(" fNECLayers %i (new) \n", fNECLayers));
 528       } else if(addOpt.Contains("PBTH=",TString::kIgnoreCase)) {//Thickness of the Pb(fECPbRadThicknes)
 529         sscanf(addOpt.Data(),"PBTH=%f", &fECPbRadThickness);
 530       } else if(addOpt.Contains("SCTH=",TString::kIgnoreCase)) {//Thickness of the Sc(fECScintThick)
 531         sscanf(addOpt.Data(),"SCTH=%f", &fECScintThick);
 532       } else if(addOpt.Contains("LATSS=",TString::kIgnoreCase)) {// Thickness of lateral steel strip (fLateralSteelStrip)
 533         sscanf(addOpt.Data(),"LATSS=%f", &fLateralSteelStrip);
 534         AliDebug(2,Form(" fLateralSteelStrip %f (new) \n", fLateralSteelStrip));
 535       }
 536     }
 537   }
 538 }
 539
 540 void AliEMCALGeometry::DefineSamplingFraction()
 541 {
 542   // Jun 05,2006
 543   // Look http://rhic.physics.wayne.edu/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
 544   // Keep for compatibilty
 545   //
 546   if(fNECLayers == 69) {        // 10% layer reduction
 547     fSampling = 12.55;
 548   } else if(fNECLayers == 61) { // 20% layer reduction
 549     fSampling = 12.80;
 550   } else if(fNECLayers == 77) {
 551     if       (fECScintThick>0.175 && fECScintThick<0.177) { // 10% Pb thicknes reduction
 552       fSampling = 10.5; // fECScintThick = 0.176, fECPbRadThickness=0.144;
 553     } else if(fECScintThick>0.191 && fECScintThick<0.193) { // 20% Pb thicknes reduction
 554       fSampling = 8.93; // fECScintThick = 0.192, fECPbRadThickness=0.128;
 555     }
 556   }
 557 }
 558
 559 //____________________________________________________________________________
 560 void AliEMCALGeometry::FillTRU(const TClonesArray * digits, TClonesArray * ampmatrix, TClonesArray * timeRmatrix) {
 561
 562
 563 //  Orders digits ampitudes list in fNTRU TRUs (384 cells) per supermodule.
 564 //  Each TRU is a TMatrixD, and they are kept in TClonesArrays. The number of
 565 //  TRU in phi is fNTRUPhi, and the number of TRU in eta is fNTRUEta.
 566 //  Last 2 modules are half size in Phi, I considered that the number of TRU
 567 //  is maintained for the last modules but decision not taken. If different,
 568 //  then this must be changed.
 569
 570
 571   //Check data members
 572
 573   if(fNTRUEta*fNTRUPhi != fNTRU)
 574     Error("FillTRU"," Wrong number of TRUS per Eta or Phi");
 575
 576   //Initilize and declare variables
 577   //List of TRU matrices initialized to 0.
 578   Int_t nCellsPhi  = fNPhi*2/fNTRUPhi;
 579   Int_t nCellsPhi2 = fNPhi/fNTRUPhi; //HalfSize modules
 580   Int_t nCellsEta  = fNZ*2/fNTRUEta;
 581   Int_t id      = -1;
 582   Float_t amp   = -1;
 583   Float_t timeR = -1;
 584   Int_t iSupMod = -1;
 585   Int_t nTower  = -1;
 586   Int_t nIphi   = -1;
 587   Int_t nIeta   = -1;
 588   Int_t iphi    = -1;
 589   Int_t ieta    = -1;
 590
 591   //List of TRU matrices initialized to 0.
 592   for(Int_t k = 0; k < fNTRU*fNumberOfSuperModules; k++){
 593     TMatrixD  * amptrus   = new TMatrixD(nCellsPhi,nCellsEta) ;
 594     TMatrixD  * timeRtrus = new TMatrixD(nCellsPhi,nCellsEta) ;
 595     for(Int_t i = 0; i < nCellsPhi; i++){
 596       for(Int_t j = 0; j < nCellsEta; j++){
 597         (*amptrus)(i,j) = 0.0;
 598         (*timeRtrus)(i,j) = 0.0;
 599       }
 600     }
 601     new((*ampmatrix)[k])   TMatrixD(*amptrus) ;
 602     new((*timeRmatrix)[k]) TMatrixD(*timeRtrus) ;
 603   }
 604
 605   AliEMCALDigit * dig ;
 606
 607   //Digits loop to fill TRU matrices with amplitudes.
 608   for(Int_t idig = 0 ; idig < digits->GetEntriesFast() ; idig++){
 609
 610     dig = dynamic_cast<AliEMCALDigit *>(digits->At(idig)) ;
 611     amp    = dig->GetAmp() ;   // Energy of the digit (arbitrary units)
 612     id     = dig->GetId() ;    // Id label of the cell
 613     timeR  = dig->GetTimeR() ; // Earliest time of the digit
 614
 615     //Get eta and phi cell position in supermodule
 616     Bool_t bCell = GetCellIndex(id, iSupMod, nTower, nIphi, nIeta) ;
 617     if(!bCell)
 618       Error("FillTRU","Wrong cell id number") ;
 619
 620     GetCellPhiEtaIndexInSModule(iSupMod,nTower,nIphi, nIeta,iphi,ieta);
 621
 622     //Check to which TRU in the supermodule belongs the cell.
 623     //Supermodules are divided in a TRU matrix of dimension
 624     //(fNTRUPhi,fNTRUEta).
 625     //Each TRU is a cell matrix of dimension (nCellsPhi,nCellsEta)
 626
 627     //First calculate the row and column in the supermodule
 628     //of the TRU to which the cell belongs.
 629     Int_t col   = ieta/nCellsEta;
 630     Int_t row   = iphi/nCellsPhi;
 631     if(iSupMod > 9)
 632       row   = iphi/nCellsPhi2;
 633     //Calculate label number of the TRU
 634     Int_t itru  = row + col*fNTRUPhi + iSupMod*fNTRU ;
 635
 636     //Fill TRU matrix with cell values
 637     TMatrixD * amptrus   = dynamic_cast<TMatrixD *>(ampmatrix->At(itru)) ;
 638     TMatrixD * timeRtrus = dynamic_cast<TMatrixD *>(timeRmatrix->At(itru)) ;
 639
 640     //Calculate row and column of the cell inside the TRU with number itru
 641     Int_t irow = iphi - row *  nCellsPhi;
 642     if(iSupMod > 9)
 643       irow = iphi - row *  nCellsPhi2;
 644     Int_t icol = ieta - col *  nCellsEta;
 645
 646     (*amptrus)(irow,icol) = amp ;
 647     (*timeRtrus)(irow,icol) = timeR ;
 648
 649   }
 650 }
 651
 652 //______________________________________________________________________
 653 void AliEMCALGeometry::GetCellPhiEtaIndexInSModuleFromTRUIndex(const Int_t itru, const Int_t iphitru, const Int_t ietatru, Int_t &iphiSM, Int_t &ietaSM) const
 654 {
 655
 656   // This method transforms the (eta,phi) index of cells in a
 657   // TRU matrix into Super Module (eta,phi) index.
 658
 659   // Calculate in which row and column where the TRU are
 660   // ordered in the SM
 661
 662   Int_t col = itru/ fNTRUPhi ;
 663   Int_t row = itru - col*fNTRUPhi ;
 664
 665   //Calculate the (eta,phi) index in SM
 666   Int_t nCellsPhi = fNPhi*2/fNTRUPhi;
 667   Int_t nCellsEta = fNZ*2/fNTRUEta;
 668
 669   iphiSM = nCellsPhi*row + iphitru  ;
 670   ietaSM = nCellsEta*col + ietatru  ;
 671 }
 672
 673 //______________________________________________________________________
 674 AliEMCALGeometry *  AliEMCALGeometry::GetInstance(){
 675   // Returns the pointer of the unique instance
 676
 677   AliEMCALGeometry * rv = static_cast<AliEMCALGeometry *>( fgGeom );
 678   return rv;
 679 }
 680
 681 //______________________________________________________________________
 682 AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
 683                                                 const Text_t* title){
 684     // Returns the pointer of the unique instance
 685
 686     AliEMCALGeometry * rv = 0;
 687     if ( fgGeom == 0 ) {
 688         if ( strcmp(name,"") == 0 ) rv = 0;
 689         else {
 690             fgGeom = new AliEMCALGeometry(name, title);
 691             if ( fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
 692             else {
 693                 rv = 0;
 694                 delete fgGeom;
 695                 fgGeom = 0;
 696             } // end if fgInit
 697         } // end if strcmp(name,"")
 698     }else{
 699         if ( strcmp(fgGeom->GetName(), name) != 0) {
 700           printf("\ncurrent geometry is %s : ", fgGeom->GetName());
 701           printf(" you cannot call %s ", name);
 702         }else{
 703           rv = (AliEMCALGeometry *) fgGeom;
 704         } // end
 705     }  // end if fgGeom
 706     return rv;
 707 }
 708
 709 Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const {
 710   // Checks whether point is inside the EMCal volume, used in AliEMCALv*.cxx
 711   //
 712   // Code uses cylindrical approximation made of inner radius (for speed)
 713   //
 714   // Points behind EMCAl, i.e. R > outer radius, but eta, phi in acceptance
 715   // are considered to inside
 716
 717   Double_t r=sqrt(x*x+y*y);
 718
 719   if ( r > fEnvelop[0] ) {
 720      Double_t theta;
 721      theta  =    TMath::ATan2(r,z);
 722      Double_t eta;
 723      if(theta == 0)
 724        eta = 9999;
 725      else
 726        eta    =   -TMath::Log(TMath::Tan(theta/2.));
 727      if (eta < fArm1EtaMin || eta > fArm1EtaMax)
 728        return 0;
 729
 730      Double_t phi = TMath::ATan2(y,x) * 180./TMath::Pi();
 731      if (phi > fArm1PhiMin && phi < fArm1PhiMax)
 732        return 1;
 733   }
 734   return 0;
 735 }
 736 // ==
 737
 738 //
 739 // == Shish-kebab cases ==
 740 //
 741 Int_t AliEMCALGeometry::GetAbsCellId(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta) const
 742 {
 743   // 27-aug-04;
 744   // corr. 21-sep-04;
 745   //       13-oct-05; 110 degree case
 746   // May 31, 2006; ALICE numbering scheme:
 747   // 0 <= nSupMod < fNumberOfSuperModules
 748   // 0 <= nTower  < fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1)
 749   // 0 <= nIphi   < fNPHIdiv
 750   // 0 <= nIeta   < fNETAdiv
 751   // 0 <= absid   < fNCells
 752   static Int_t id=0; // have to change from 0 to fNCells-1
 753   if(fKey110DEG == 1 && nSupMod >= 10) { // 110 degree case; last two supermodules
 754     id  = fNCellsInSupMod*10 + (fNCellsInSupMod/2)*(nSupMod-10);
 755   } else {
 756     id  = fNCellsInSupMod*nSupMod;
 757   }
 758   id += fNCellsInTower *nTower;
 759   id += fNPHIdiv *nIphi;
 760   id += nIeta;
 761   if(id<0 || id >= fNCells) {
 762 //     printf(" wrong numerations !!\n");
 763 //     printf("    id      %6i(will be force to -1)\n", id);
 764 //     printf("    fNCells %6i\n", fNCells);
 765 //     printf("    nSupMod %6i\n", nSupMod);
 766 //     printf("    nTower  %6i\n", nTower);
 767 //     printf("    nIphi   %6i\n", nIphi);
 768 //     printf("    nIeta   %6i\n", nIeta);
 769     id = -TMath::Abs(id); // if negative something wrong
 770   }
 771   return id;
 772 }
 773
 774 Bool_t  AliEMCALGeometry::CheckAbsCellId(Int_t absId) const
 775 {
 776   // May 31, 2006; only trd1 now
 777   if(absId<0 || absId >= fNCells) return kFALSE;
 778   else                            return kTRUE;
 779 }
 780
 781 Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nTower,Int_t &nIphi,Int_t &nIeta) const
 782 {
 783   // 21-sep-04; 19-oct-05;
 784   // May 31, 2006; ALICE numbering scheme:
 785   static Int_t tmp=0, sm10=0;
 786   if(!CheckAbsCellId(absId)) return kFALSE;
 787
 788   sm10 = fNCellsInSupMod*10;
 789   if(fKey110DEG == 1 && absId >= sm10) { // 110 degree case; last two supermodules
 790     nSupMod = (absId-sm10) / (fNCellsInSupMod/2) + 10;
 791     tmp     = (absId-sm10) % (fNCellsInSupMod/2);
 792   } else {
 793     nSupMod = absId / fNCellsInSupMod;
 794     tmp     = absId % fNCellsInSupMod;
 795   }
 796
 797   nTower  = tmp / fNCellsInTower;
 798   tmp     = tmp % fNCellsInTower;
 799   nIphi   = tmp / fNPHIdiv;
 800   nIeta   = tmp % fNPHIdiv;
 801
 802   return kTRUE;
 803 }
 804
 805 void AliEMCALGeometry::GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower,  int &iphim, int &ietam) const
 806 {
 807   // added nSupMod; - 19-oct-05 !
 808   // Alice numbering scheme        - Jun 01,2006
 809   // ietam, iphi - indexes of module in two dimensional grid of SM
 810   // ietam - have to change from 0 to fNZ-1
 811   // iphim - have to change from 0 to nphi-1 (fNPhi-1 or fNPhi/2-1)
 812   static Int_t nphi;
 813
 814   if(fKey110DEG == 1 && nSupMod>=10) nphi = fNPhi/2;
 815   else                               nphi = fNPhi;
 816
 817   ietam = nTower/nphi;
 818   iphim = nTower%nphi;
 819 }
 820
 821 void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta,
 822 int &iphi, int &ieta) const
 823 {
 824   //
 825   // Added nSupMod; Nov 25, 05
 826   // Alice numbering scheme  - Jun 01,2006
 827   // ieta, iphi - indexes of cell(tower) in two dimensional grid of SM
 828   // ieta - have to change from 0 to (fNZ*fNETAdiv-1)
 829   // iphi - have to change from 0 to (fNPhi*fNPHIdiv-1 or fNPhi*fNPHIdiv/2-1)
 830   //
 831   static Int_t iphim, ietam;
 832
 833   GetModulePhiEtaIndexInSModule(nSupMod,nTower, iphim, ietam);
 834   //  ieta  = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM)
 835   ieta  = ietam*fNETAdiv + (fNETAdiv - 1 - nIeta); // x(module) = -z(SM)
 836   iphi  = iphim*fNPHIdiv + nIphi;     // y(module) =  y(SM)
 837
 838   if(iphi<0 || ieta<0)
 839   AliDebug(1,Form(" nSupMod %i nTower %i nIphi %i nIeta %i => ieta %i iphi %i\n",
 840   nSupMod, nTower, nIphi, nIeta, ieta, iphi));
 841 }
 842
 843 Int_t  AliEMCALGeometry::GetSuperModuleNumber(Int_t absId)  const
 844 {
 845   // Return the number of the  supermodule given the absolute
 846   // ALICE numbering id
 847
 848   static Int_t nSupMod, nTower, nIphi, nIeta;
 849   GetCellIndex(absId, nSupMod, nTower, nIphi, nIeta);
 850   return nSupMod;
 851 }
 852
 853 void  AliEMCALGeometry::GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta,
 854                         Int_t &iphim, Int_t &ietam, Int_t &nTower) const
 855 {
 856   // Transition from cell indexes (ieta,iphi) to module indexes (ietam,iphim, nTower)
 857   static Int_t nphi;
 858   nphi  = GetNumberOfModuleInPhiDirection(nSupMod);
 859
 860   ietam  = ieta/fNETAdiv;
 861   iphim  = iphi/fNPHIdiv;
 862   nTower = ietam * nphi + iphim;
 863 }
 864
 865 Int_t  AliEMCALGeometry::GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const
 866 {
 867   // Transition from super module number(nSupMod) and cell indexes (ieta,iphi) to absId
 868   static Int_t ietam, iphim, nTower;
 869   static Int_t nIeta, nIphi; // cell indexes in module
 870
 871   GetModuleIndexesFromCellIndexesInSModule(nSupMod, iphi, ieta, ietam, iphim, nTower);
 872
 873   nIeta = ieta%fNETAdiv;
 874   nIeta = fNETAdiv - 1 - nIeta;
 875   nIphi = iphi%fNPHIdiv;
 876
 877   return GetAbsCellId(nSupMod, nTower, nIphi, nIeta);
 878 }
 879
 880
 881 // Methods for AliEMCALRecPoint - Feb 19, 2006
 882 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const
 883 {
 884   // Look to see what the relative
 885   // position inside a given cell is
 886   // for a recpoint.
 887   // Alice numbering scheme - Jun 08, 2006
 888
 889   // Shift index taking into account the difference between standard SM
 890   // and SM of half size in phi direction
 891   const Int_t phiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
 892   static Int_t nSupMod, nTower, nIphi, nIeta, iphi, ieta;
 893   if(!CheckAbsCellId(absId)) return kFALSE;
 894
 895   GetCellIndex(absId, nSupMod, nTower, nIphi, nIeta);
 896   GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi, ieta);
 897
 898   xr = fCentersOfCellsXDir.At(ieta);
 899   zr = fCentersOfCellsEtaDir.At(ieta);
 900
 901   if(nSupMod<10) {
 902     yr = fCentersOfCellsPhiDir.At(iphi);
 903   } else {
 904     yr = fCentersOfCellsPhiDir.At(iphi + phiIndexShift);
 905   }
 906   AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
 907
 908   return kTRUE;
 909 }
 910
 911 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t loc[3]) const
 912 {
 913   // Alice numbering scheme - Jun 03, 2006
 914   loc[0] = loc[1] = loc[2]=0.0;
 915   if(RelPosCellInSModule(absId, loc[0],loc[1],loc[2])) {
 916     return kTRUE;
 917   }
 918   return kFALSE;
 919 }
 920
 921 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, TVector3 &vloc) const
 922 {
 923   static Double_t loc[3];
 924   if(RelPosCellInSModule(absId,loc)) {
 925     vloc.SetXYZ(loc[0], loc[1], loc[2]);
 926     return kTRUE;
 927   } else {
 928     vloc.SetXYZ(0,0,0);
 929     return kFALSE;
 930   }
 931   // Alice numbering scheme - Jun 03, 2006
 932 }
 933
 934 void AliEMCALGeometry::CreateListOfTrd1Modules()
 935 {
 936   // Generate the list of Trd1 modules
 937   // which will make up the EMCAL
 938   // geometry
 939
 940   AliDebug(2,Form(" AliEMCALGeometry::CreateListOfTrd1Modules() started "));
 941
 942   AliEMCALShishKebabTrd1Module *mod=0, *mTmp=0; // current module
 943   if(fShishKebabTrd1Modules == 0) {
 944     fShishKebabTrd1Modules = new TList;
 945     fShishKebabTrd1Modules->SetName("ListOfTRD1");
 946     for(int iz=0; iz< GetNZ(); iz++) {
 947       if(iz==0) {
 948         mod  = new AliEMCALShishKebabTrd1Module(TMath::Pi()/2.,this);
 949       } else {
 950         mTmp  = new AliEMCALShishKebabTrd1Module(*mod);
 951         mod   = mTmp;
 952       }
 953       fShishKebabTrd1Modules->Add(mod);
 954     }
 955   } else {
 956     AliDebug(2,Form(" Already exits : "));
 957   }
 958   mod = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(fShishKebabTrd1Modules->GetSize()-1);
 959   fEtaMaxOfTRD1 = mod->GetMaxEtaOfModule(0);
 960
 961   AliDebug(2,Form(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n",
 962                   fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1));
 963   // Feb 20,2006;
 964   // Jun 01, 2006 - ALICE numbering scheme
 965   // define grid for cells in eta(z) and x directions in local coordinates system of SM
 966   // Works just for 2x2 case only -- ?? start here
 967   //
 968   //
 969   // Define grid for cells in phi(y) direction in local coordinates system of SM
 970   // as for 2X2 as for 3X3 - Nov 8,2006
 971   //
 972   AliDebug(2,Form(" Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize()));
 973   Int_t ind=0; // this is phi index
 974   Int_t iphi=0, ieta=0, nTower=0, iphiTemp;
 975   Double_t xr, zr, theta, phi, eta, r, x,y;
 976   TVector3 vglob;
 977   Double_t ytCenterModule, ytCenterCell;
 978
 979   fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
 980   fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
 981
 982   Double_t R0 = GetIPDistance() + GetLongModuleSize()/2.;
 983   for(Int_t it=0; it<fNPhi; it++) { // cycle on modules
 984     ytCenterModule = -fParSM[1] + fPhiModuleSize*(2*it+1)/2;  // center of module
 985     for(Int_t ic=0; ic<fNPHIdiv; ic++) { // cycle on cells in module
 986       if(fNPHIdiv==2) {
 987         ytCenterCell = ytCenterModule + fPhiTileSize *(2*ic-1)/2.;
 988       } else if(fNPHIdiv==3){
 989         ytCenterCell = ytCenterModule + fPhiTileSize *(ic-1);
 990       } else if(fNPHIdiv==1){
 991         ytCenterCell = ytCenterModule;
 992       }
 993       fCentersOfCellsPhiDir.AddAt(ytCenterCell,ind);
 994       // Define grid on phi direction
 995       // Grid is not the same for different eta bin;
 996       // Effect is small but is still here
 997       phi = TMath::ATan2(ytCenterCell, R0);
 998       fPhiCentersOfCells.AddAt(phi, ind);
 999
1000       AliDebug(2,Form(" ind %2.2i : y %8.3f ", ind, fCentersOfCellsPhiDir.At(ind)));
1001       ind++;
1002     }
1003   }
1004
1005   fCentersOfCellsEtaDir.Set(fNZ *fNETAdiv);
1006   fCentersOfCellsXDir.Set(fNZ *fNETAdiv);
1007   fEtaCentersOfCells.Set(fNZ *fNETAdiv * fNPhi*fNPHIdiv);
1008   AliDebug(2,Form(" Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize()));
1009   for(Int_t it=0; it<fNZ; it++) {
1010     AliEMCALShishKebabTrd1Module *trd1 = GetShishKebabModule(it);
1011     nTower = fNPhi*it;
1012     for(Int_t ic=0; ic<fNETAdiv; ic++) {
1013       if(fNPHIdiv==2) {
1014         trd1->GetCenterOfCellInLocalCoordinateofSM(ic, xr, zr);      // case of 2X2
1015         GetCellPhiEtaIndexInSModule(0, nTower, 0, ic, iphiTemp, ieta);
1016       } if(fNPHIdiv==3) {
1017         trd1->GetCenterOfCellInLocalCoordinateofSM_3X3(ic, xr, zr);  // case of 3X3
1018         GetCellPhiEtaIndexInSModule(0, nTower, 0, ic, iphiTemp, ieta);
1019       } if(fNPHIdiv==1) {
1020         trd1->GetCenterOfCellInLocalCoordinateofSM_1X1(xr, zr);      // case of 1X1
1021         GetCellPhiEtaIndexInSModule(0, nTower, 0, ic, iphiTemp, ieta);
1022       }
1023       fCentersOfCellsXDir.AddAt(float(xr) - fParSM[0],ieta);
1024       fCentersOfCellsEtaDir.AddAt(float(zr) - fParSM[2],ieta);
1025       // Define grid on eta direction for each bin in phi
1026       for(int iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
1027         x = xr + trd1->GetRadius();
1028         y = fCentersOfCellsPhiDir[iphi];
1029         r = TMath::Sqrt(x*x + y*y + zr*zr);
1030         theta = TMath::ACos(zr/r);
1031         eta   = AliEMCALShishKebabTrd1Module::ThetaToEta(theta);
1032         //        ind   = ieta*fCentersOfCellsPhiDir.GetSize() + iphi;
1033         ind   = iphi*fCentersOfCellsEtaDir.GetSize() + ieta;
1034         fEtaCentersOfCells.AddAt(eta, ind);
1035       }
1036       //printf(" ieta %i : xr + trd1->GetRadius() %f : zr %f : eta %f \n", ieta, xr + trd1->GetRadius(), zr, eta);
1037     }
1038   }
1039   for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1040     AliDebug(2,Form(" ind %2.2i : z %8.3f : x %8.3f", i+1,
1041                     fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i)));
1042   }
1043
1044 }
1045
1046 void  AliEMCALGeometry::GetTransformationForSM()
1047 {
1048   //Uses the geometry manager to
1049   //load the transformation matrix
1050   //for the supermodules
1051
1052   static Bool_t transInit=kFALSE;
1053   if(transInit) return;
1054
1055   int i=0;
1056   if(gGeoManager == 0) {
1057     Info("CreateTransformationForSM() "," Load geometry : TGeoManager::Import()");
1058     assert(0);
1059   }
1060   TGeoNode *tn = gGeoManager->GetTopNode();
1061   TGeoNode *node=0, *xen1 = 0;
1062   for(i=0; i<tn->GetNdaughters(); i++) {
1063     node = tn->GetDaughter(i);
1064     TString ns(node->GetName());
1065     if(ns.Contains(GetNameOfEMCALEnvelope())) {
1066       xen1 = node;
1067       break;
1068     }
1069   }
1070   if(!xen1) {
1071     Info("CreateTransformationForSM() "," geometry has not EMCAL envelope with name %s",
1072     GetNameOfEMCALEnvelope());
1073     assert(0);
1074   }
1075   printf(" i %i : EMCAL Envelope is %s : #SM %i \n", i, xen1->GetName(), xen1->GetNdaughters());
1076   for(i=0; i<xen1->GetNdaughters(); i++) {
1077     TGeoNodeMatrix *sm = (TGeoNodeMatrix*)xen1->GetDaughter(i);
1078     fMatrixOfSM[i] = sm->GetMatrix();
1079     //Compiler doesn't like this syntax...
1080     //    printf(" %i : matrix %x \n", i, fMatrixOfSM[i]);
1081   }
1082   transInit = kTRUE;
1083 }
1084
1085 void AliEMCALGeometry::GetGlobal(const Double_t *loc, Double_t *glob, int ind) const
1086 {
1087   // Figure out the global numbering
1088   // of a given supermodule from the
1089   // local numbering
1090   // Alice numbering - Jun 03,2006
1091   //  if(fMatrixOfSM[0] == 0) GetTransformationForSM();
1092
1093   if(ind>=0 && ind < GetNumberOfSuperModules()) {
1094     fMatrixOfSM[ind]->LocalToMaster(loc, glob);
1095   }
1096 }
1097
1098 void AliEMCALGeometry::GetGlobal(const TVector3 &vloc, TVector3 &vglob, int ind) const
1099 {
1100   //Figure out the global numbering
1101   //of a given supermodule from the
1102   //local numbering given a 3-vector location
1103
1104   static Double_t tglob[3], tloc[3];
1105   vloc.GetXYZ(tloc);
1106   GetGlobal(tloc, tglob, ind);
1107   vglob.SetXYZ(tglob[0], tglob[1], tglob[2]);
1108 }
1109
1110 void AliEMCALGeometry::GetGlobal(Int_t absId , double glob[3]) const
1111 {
1112   // Alice numbering scheme - Jun 03, 2006
1113   static Int_t nSupMod, nModule, nIphi, nIeta;
1114   static double loc[3];
1115
1116   glob[0]=glob[1]=glob[2]=0.0; // bad case
1117   if(RelPosCellInSModule(absId, loc)) {
1118     GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1119     fMatrixOfSM[nSupMod]->LocalToMaster(loc, glob);
1120   }
1121 }
1122
1123 void AliEMCALGeometry::GetGlobal(Int_t absId , TVector3 &vglob) const
1124 {
1125   // Alice numbering scheme - Jun 03, 2006
1126   static Double_t glob[3];
1127
1128   GetGlobal(absId, glob);
1129   vglob.SetXYZ(glob[0], glob[1], glob[2]);
1130
1131 }
1132
1133 void AliEMCALGeometry::GetGlobal(const AliRecPoint *rp, TVector3 &vglob) const
1134 {
1135   // Figure out the global numbering
1136   // of a given supermodule from the
1137   // local numbering for RecPoints
1138
1139   static TVector3 vloc;
1140   static Int_t nSupMod, nModule, nIphi, nIeta;
1141
1142   AliRecPoint *rpTmp = (AliRecPoint*)rp; // const_cast ??
1143   if(!rpTmp) return;
1144   AliEMCALRecPoint *rpEmc = (AliEMCALRecPoint*)rpTmp;
1145
1146   GetCellIndex(rpEmc->GetAbsId(0), nSupMod, nModule, nIphi, nIeta);
1147   rpTmp->GetLocalPosition(vloc);
1148   GetGlobal(vloc, vglob, nSupMod);
1149 }
1150
1151 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Double_t &eta,Double_t &phi) const
1152 {
1153   // Nov 16, 2006- float to double
1154   // version for TRD1 only
1155   static TVector3 vglob;
1156   GetGlobal(absId, vglob);
1157   eta = vglob.Eta();
1158   phi = vglob.Phi();
1159 }
1160
1161 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Float_t &eta,Float_t &phi) const
1162 {
1163   // Nov 16,2006 - should be discard in future
1164   static TVector3 vglob;
1165   GetGlobal(absId, vglob);
1166   eta = float(vglob.Eta());
1167   phi = float(vglob.Phi());
1168 }
1169
1170 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSM(Int_t nSupMod, Double_t &phiMin, Double_t &phiMax) const
1171 {
1172   // 0<= nSupMod <=11; phi in rad
1173   static int i;
1174   if(nSupMod<0 || nSupMod >11) return kFALSE;
1175   i = nSupMod/2;
1176   phiMin = fPhiBoundariesOfSM[2*i];
1177   phiMax = fPhiBoundariesOfSM[2*i+1];
1178   return kTRUE;
1179 }
1180
1181 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const
1182 {
1183   // 0<= nPhiSec <=4; phi in rad
1184   // 0;  gap boundaries between  0th&2th  | 1th&3th SM
1185   // 1;  gap boundaries between  2th&4th  | 3th&5th SM
1186   // 2;  gap boundaries between  4th&6th  | 5th&7th SM
1187   // 3;  gap boundaries between  6th&8th  | 7th&9th SM
1188   // 4;  gap boundaries between  8th&10th | 9th&11th SM
1189   if(nPhiSec<0 || nPhiSec >4) return kFALSE;
1190   phiMin = fPhiBoundariesOfSM[2*nPhiSec+1];
1191   phiMax = fPhiBoundariesOfSM[2*nPhiSec+2];
1192   return kTRUE;
1193 }
1194
1195 Bool_t AliEMCALGeometry::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const
1196 {
1197   // Return false if phi belongs a phi cracks between SM
1198
1199   static Int_t i;
1200
1201   if(TMath::Abs(eta) > fEtaMaxOfTRD1) return kFALSE;
1202
1203   phi = TVector2::Phi_0_2pi(phi); // move phi to (0,2pi) boundaries
1204   for(i=0; i<6; i++) {
1205     if(phi>=fPhiBoundariesOfSM[2*i] && phi<=fPhiBoundariesOfSM[2*i+1]) {
1206       nSupMod = 2*i;
1207       if(eta < 0.0) nSupMod++;
1208       AliDebug(1,Form("eta %f phi %f(%5.2f) : nSupMod %i : #bound %i", eta,phi,phi*TMath::RadToDeg(), nSupMod,i));
1209       return kTRUE;
1210     }
1211   }
1212   return kFALSE;
1213 }
1214
1215 Bool_t AliEMCALGeometry::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi, Int_t &absId) const
1216 {
1217   // Nov 17,2006
1218   // stay here - phi problem as usual
1219   static Int_t nSupMod, i, ieta, iphi, etaShift, nphi;
1220   static Double_t absEta=0.0, d=0.0, dmin=0.0, phiLoc;
1221   absId = nSupMod = - 1;
1222   if(SuperModuleNumberFromEtaPhi(eta, phi, nSupMod)) {
1223     // phi index first
1224     phi    = TVector2::Phi_0_2pi(phi);
1225     phiLoc = phi - fPhiCentersOfSM[nSupMod/2];
1226     nphi   = fPhiCentersOfCells.GetSize();
1227     if(nSupMod>=10) {
1228       phiLoc = phi - 190.*TMath::DegToRad();
1229       nphi  /= 2;
1230     }
1231
1232     dmin   = TMath::Abs(fPhiCentersOfCells[0]-phiLoc);
1233     iphi   = 0;
1234     for(i=1; i<nphi; i++) {
1235       d = TMath::Abs(fPhiCentersOfCells[i] - phiLoc);
1236       if(d < dmin) {
1237         dmin = d;
1238         iphi = i;
1239       }
1240       //      printf(" i %i : d %f : dmin %f : fPhiCentersOfCells[i] %f \n", i, d, dmin, fPhiCentersOfCells[i]);
1241     }
1242     // odd SM are turned with respect of even SM - reverse indexes
1243     AliDebug(2,Form(" iphi %i : dmin %f (phi %f, phiLoc %f ) ", iphi, dmin, phi, phiLoc));
1244     // eta index
1245     absEta   = TMath::Abs(eta);
1246     etaShift = iphi*fCentersOfCellsEtaDir.GetSize();
1247     dmin     = TMath::Abs(fEtaCentersOfCells[etaShift]-absEta);
1248     ieta     = 0;
1249     for(i=1; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1250       d = TMath::Abs(fEtaCentersOfCells[i+etaShift] - absEta);
1251       if(d < dmin) {
1252         dmin = d;
1253         ieta = i;
1254       }
1255     }
1256     AliDebug(2,Form(" ieta %i : dmin %f (eta=%f) : nSupMod %i ", ieta, dmin, eta, nSupMod));
1257
1258     if(eta<0) iphi = (nphi-1) - iphi;
1259     absId = GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);
1260
1261     return kTRUE;
1262   }
1263   return kFALSE;
1264 }
1265
1266 AliEMCALShishKebabTrd1Module* AliEMCALGeometry::GetShishKebabModule(Int_t neta)
1267 {
1268   //This method was too long to be
1269   //included in the header file - the
1270   //rule checker complained about it's
1271   //length, so we move it here.  It returns the
1272   //shishkebabmodule at a given eta index point.
1273
1274   static AliEMCALShishKebabTrd1Module* trd1=0;
1275   if(fShishKebabTrd1Modules && neta>=0 && neta<fShishKebabTrd1Modules->GetSize()) {
1276     trd1 = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(neta);
1277   } else trd1 = 0;
1278   return trd1;
1279 }
1280
1281 void AliEMCALGeometry::Browse(TBrowser* b)
1282 {
1283   if(fShishKebabTrd1Modules) b->Add(fShishKebabTrd1Modules);
1284 }
1285
1286 Bool_t AliEMCALGeometry::IsFolder() const
1287 {
1288   if(fShishKebabTrd1Modules) return kTRUE;
1289   else                       return kFALSE;
1290 }