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