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