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