EMCAL/AliEMCALGeometry.cxx

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