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