EMCAL/AliEMCALGeometry.cxx

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