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