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