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 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::GetRunLoader();
  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"))) {
 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   // constant for transition absid <--> indexes
 335   fNCellsInModule  = fNPHIdiv*fNETAdiv;
 336   fNCellsInSupMod = fNCellsInModule*fNPhi*fNZ;
 337   fNCells         = fNCellsInSupMod*fNumberOfSuperModules;
 338   if(GetKey110DEG()) fNCells -= fNCellsInSupMod;
 339
 340   fNPhiSuperModule = fNumberOfSuperModules/2;
 341   if(fNPhiSuperModule<1) fNPhiSuperModule = 1;
 342
 343   fPhiTileSize = fPhiModuleSize/double(fNPHIdiv) - fLateralSteelStrip; // 13-may-05
 344   fEtaTileSize = fEtaModuleSize/double(fNETAdiv) - fLateralSteelStrip; // 13-may-05
 345
 346   fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness);
 347   f2Trd1Dx2 = fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.);
 348   if(!fGeoName.Contains("WSUC")) fShellThickness  = TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2);
 349
 350   //These parameters are used to create the mother volume to hold the supermodules
 351   //2cm padding added to allow for misalignments - JLK 30-May-2008
 352   fEnvelop[0]     = fIPDistance - 1.; // mother volume inner radius
 353   fEnvelop[1]     = fIPDistance + fShellThickness + 1.; // mother volume outer r.
 354   fEnvelop[2]     = fZLength + 2.; //mother volume length
 355
 356   // Local coordinates
 357   fParSM[0] = GetShellThickness()/2.;
 358   fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
 359   fParSM[2] = fZLength/4.;  //divide by 4 to get half-length of SM
 360
 361   // SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006
 362   fPhiBoundariesOfSM.Set(fNumberOfSuperModules);
 363   fPhiCentersOfSM.Set(fNumberOfSuperModules/2);
 364   fPhiBoundariesOfSM[0] = TMath::PiOver2() - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules)
 365   fPhiCentersOfSM[0]     = TMath::PiOver2();
 366   if(fNumberOfSuperModules > 1)
 367     fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance);
 368   if(fNumberOfSuperModules > 2) {
 369     for(int i=1; i<=4; i++) { // from 2th ro 9th
 370       fPhiBoundariesOfSM[2*i]   = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i;
 371       fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i;
 372       fPhiCentersOfSM[i]         = fPhiCentersOfSM[0]     + 20.*TMath::DegToRad()*i;
 373     }
 374   }
 375   if(fNumberOfSuperModules > 10) {
 376     fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad();
 377     fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance);
 378     fPhiCentersOfSM[5]      = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.;
 379   }
 380
 381   //called after setting of scintillator and lead layer parameters
 382   DefineSamplingFraction();
 383
 384   // TRU parameters - Apr 29,08 by PAI.
 385   // These parameters values was updated at Nov 05, 2007
 386   // As is on Olivier  BOURRION (LPSC) ppt preasentation
 387   // at ALICE trigger meeting at 13th-14th March
 388   fNTRUEta = 1;           // was 3
 389   fNTRUPhi = 3;           // was 1
 390   fNModulesInTRUEta = 24; // was 8
 391   fNModulesInTRUPhi = 4;  // was 12
 392   // Jet trigger
 393   // 3*6*10 + 2*6*2 = 204 -> matrix (nphi(17), neta(12))
 394   fNEtaSubOfTRU     = 6;
 395
 396   fgInit = kTRUE;
 397 }
 398
 399 void AliEMCALGeometry::PrintGeometry()
 400 {
 401   // Separate routine is callable from broswer; Nov 7,2006
 402   printf("\nInit: geometry of EMCAL named %s :\n", fGeoName.Data());
 403   if(fArrayOpts) {
 404     for(Int_t i=0; i<fArrayOpts->GetEntries(); i++){
 405       TObjString *o = (TObjString*)fArrayOpts->At(i);
 406       printf(" %i : %s \n", i, o->String().Data());
 407     }
 408   }
 409   printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
 410   printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f -> for EMCAL envelope only\n",
 411            GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() );
 412
 413   printf( "               ECAL      : %d x (%f cm Pb, %f cm Sc) \n",
 414   GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ;
 415   printf("                fSampling %5.2f \n",  fSampling );
 416   printf(" fIPDistance       %6.3f cm \n", fIPDistance);
 417   printf(" fNPhi %i   |  fNZ %i \n", fNPhi, fNZ);
 418   printf(" fNCellsInModule %i : fNCellsInSupMod %i : fNCells %i\n",fNCellsInModule, fNCellsInSupMod, fNCells);
 419   printf(" X:Y module size     %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize);
 420   printf(" X:Y   tile size     %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize);
 421   printf(" #of sampling layers %i(fNECLayers) \n", fNECLayers);
 422   printf(" fLongModuleSize     %6.3f cm \n", fLongModuleSize);
 423   printf(" #supermodule in phi direction %i \n", fNPhiSuperModule );
 424   printf(" fILOSS %i : fIHADR %i \n", fILOSS, fIHADR);
 425   printf(" fTrd1Angle %7.4f\n", fTrd1Angle);
 426   printf(" f2Trd1Dx2  %7.4f\n",  f2Trd1Dx2);
 427   printf("SM dimensions(TRD1) : dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n",
 428          fParSM[0],fParSM[1],fParSM[2]);
 429   printf(" fPhiGapForSM  %7.4f cm (%7.4f <- phi size in degree)\n",
 430          fPhiGapForSM, TMath::ATan2(fPhiGapForSM,fIPDistance)*TMath::RadToDeg());
 431   if(GetKey110DEG()) printf(" Last two modules have size 10 degree in  phi (180<phi<190)\n");
 432   printf(" phi SM boundaries \n");
 433   for(int i=0; i<fPhiBoundariesOfSM.GetSize()/2.; i++) {
 434     printf(" %i : %7.5f(%7.2f) -> %7.5f(%7.2f) : center %7.5f(%7.2f) \n", i,
 435            fPhiBoundariesOfSM[2*i], fPhiBoundariesOfSM[2*i]*TMath::RadToDeg(),
 436            fPhiBoundariesOfSM[2*i+1], fPhiBoundariesOfSM[2*i+1]*TMath::RadToDeg(),
 437            fPhiCentersOfSM[i], fPhiCentersOfSM[i]*TMath::RadToDeg());
 438   }
 439   printf(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n",
 440          fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1);
 441
 442   printf("\n Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize());
 443   for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
 444     printf(" ind %2.2i : z %8.3f : x %8.3f \n", i,
 445            fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i));
 446     int ind=0; // Nov 21,2006
 447     for(Int_t iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
 448       ind = iphi*fCentersOfCellsEtaDir.GetSize() + i;
 449       printf("%6.4f ", fEtaCentersOfCells[ind]);
 450       if((iphi+1)%12 == 0) printf("\n");
 451     }
 452     printf("\n");
 453
 454   }
 455
 456   printf("\n Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize());
 457   for(Int_t i=0; i<fCentersOfCellsPhiDir.GetSize(); i++) {
 458     double phi=fPhiCentersOfCells.At(i);
 459     printf(" ind %2.2i : y %8.3f : phi %7.5f(%6.2f) \n", i, fCentersOfCellsPhiDir.At(i),
 460            phi, phi*TMath::RadToDeg());
 461   }
 462
 463 }
 464
 465 //______________________________________________________________________
 466 void AliEMCALGeometry::PrintCellIndexes(Int_t absId, int pri, char *tit)
 467 {
 468   // Service methods
 469   Int_t nSupMod, nModule, nIphi, nIeta;
 470   Int_t iphi, ieta;
 471   TVector3 vg;
 472
 473   GetCellIndex(absId,  nSupMod, nModule, nIphi, nIeta);
 474   printf(" %s | absId : %i -> nSupMod %i nModule %i nIphi %i nIeta %i \n", tit, absId,  nSupMod, nModule, nIphi, nIeta);
 475   if(pri>0) {
 476     GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
 477     printf(" local SM index : iphi %i : ieta %i \n", iphi,ieta);
 478     GetGlobal(absId, vg);
 479     printf(" vglob : mag %7.2f : perp %7.2f : z %7.2f : eta %6.4f : phi %6.4f(%6.2f) \n",
 480            vg.Mag(), vg.Perp(), vg.Z(), vg.Eta(), vg.Phi(), vg.Phi()*TMath::RadToDeg());
 481   }
 482 }
 483
 484 //______________________________________________________________________
 485 void AliEMCALGeometry::CheckAdditionalOptions()
 486 {
 487   // Feb 06,2006
 488   // Additional options that
 489   // can be used to select
 490   // the specific geometry of
 491   // EMCAL to run
 492   // Dec 27,2006
 493   // adeed allILOSS= and allIHADR= for MIP investigation
 494   fArrayOpts = new TObjArray;
 495   Int_t nopt = AliEMCALHistoUtilities::ParseString(fGeoName, *fArrayOpts);
 496   if(nopt==1) { // no aditional option(s)
 497     fArrayOpts->Delete();
 498     delete fArrayOpts;
 499     fArrayOpts = 0;
 500     return;
 501   }
 502   for(Int_t i=1; i<nopt; i++){
 503     TObjString *o = (TObjString*)fArrayOpts->At(i);
 504
 505     TString addOpt = o->String();
 506     Int_t indj=-1;
 507     for(Int_t j=0; j<fNAdditionalOpts; j++) {
 508       TString opt = fAdditionalOpts[j];
 509       if(addOpt.Contains(opt,TString::kIgnoreCase)) {
 510           indj = j;
 511         break;
 512       }
 513     }
 514     if(indj<0) {
 515       AliDebug(2,Form("<E> option |%s| unavailable : ** look to the file AliEMCALGeometry.h **\n",
 516                       addOpt.Data()));
 517       assert(0);
 518     } else {
 519       AliDebug(2,Form("<I> option |%s| is valid : number %i : |%s|\n",
 520                       addOpt.Data(), indj, fAdditionalOpts[indj]));
 521       if       (addOpt.Contains("NL=",TString::kIgnoreCase))   {// number of sampling layers
 522         sscanf(addOpt.Data(),"NL=%i", &fNECLayers);
 523         AliDebug(2,Form(" fNECLayers %i (new) \n", fNECLayers));
 524       } else if(addOpt.Contains("PBTH=",TString::kIgnoreCase)) {//Thickness of the Pb(fECPbRadThicknes)
 525         sscanf(addOpt.Data(),"PBTH=%f", &fECPbRadThickness);
 526       } else if(addOpt.Contains("SCTH=",TString::kIgnoreCase)) {//Thickness of the Sc(fECScintThick)
 527         sscanf(addOpt.Data(),"SCTH=%f", &fECScintThick);
 528       } else if(addOpt.Contains("LATSS=",TString::kIgnoreCase)) {// Thickness of lateral steel strip (fLateralSteelStrip)
 529         sscanf(addOpt.Data(),"LATSS=%f", &fLateralSteelStrip);
 530         AliDebug(2,Form(" fLateralSteelStrip %f (new) \n", fLateralSteelStrip));
 531       } else if(addOpt.Contains("ILOSS=",TString::kIgnoreCase)) {// As in Geant
 532         sscanf(addOpt.Data(),"ALLILOSS=%i", &fILOSS);
 533         AliDebug(2,Form(" fILOSS %i \n", fILOSS));
 534       } else if(addOpt.Contains("IHADR=",TString::kIgnoreCase)) {// As in Geant
 535         sscanf(addOpt.Data(),"ALLIHADR=%i", &fIHADR);
 536         AliDebug(2,Form(" fIHADR %i \n", fIHADR));
 537       }
 538     }
 539   }
 540 }
 541
 542 void AliEMCALGeometry::DefineSamplingFraction()
 543 {
 544   // Jun 05,2006
 545   // Look http://rhic.physics.wayne.edu/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
 546   // Keep for compatibilty
 547   //
 548   if(fNECLayers == 69) {        // 10% layer reduction
 549     fSampling = 12.55;
 550   } else if(fNECLayers == 61) { // 20% layer reduction
 551     fSampling = 12.80;
 552   } else if(fNECLayers == 77) {
 553     if       (fECScintThick>0.159 && fECScintThick<0.161) { // original sampling fraction, equal layers
 554       fSampling = 12.327; // fECScintThick = fECPbRadThickness = 0.160;
 555     } else if (fECScintThick>0.175 && fECScintThick<0.177) { // 10% Pb thicknes reduction
 556       fSampling = 10.5; // fECScintThick = 0.176, fECPbRadThickness=0.144;
 557     } else if(fECScintThick>0.191 && fECScintThick<0.193) { // 20% Pb thicknes reduction
 558       fSampling = 8.93; // fECScintThick = 0.192, fECPbRadThickness=0.128;
 559     }
 560
 561   }
 562 }
 563
 564 //______________________________________________________________________
 565 void AliEMCALGeometry::GetModulePhiEtaIndexInSModuleFromTRUIndex(const Int_t itru, const Int_t iphitru, const Int_t ietatru, Int_t &iphiSM, Int_t &ietaSM) const
 566 {
 567
 568   // This method transforms the (eta,phi) index of module in a
 569   // TRU matrix into Super Module (eta,phi) index.
 570
 571   // Calculate in which row and column where the TRU are
 572   // ordered in the SM
 573
 574   Int_t col = itru/ fNTRUPhi ; // indexes of TRU in SM
 575   Int_t row = itru - col*fNTRUPhi ;
 576
 577   iphiSM = fNModulesInTRUPhi*row + iphitru  ;
 578   ietaSM = fNModulesInTRUEta*col + ietatru  ;
 579   //printf(" GetModulePhiEtaIndexInSModuleFromTRUIndex : itru %2i iphitru %2i ietatru %2i iphiSM %2i ietaSM %2i \n",
 580   // itru, iphitru, ietatru, iphiSM, ietaSM);
 581 }
 582
 583 //______________________________________________________________________
 584 AliEMCALGeometry *  AliEMCALGeometry::GetInstance(){
 585   // Returns the pointer of the unique instance
 586
 587   AliEMCALGeometry * rv = static_cast<AliEMCALGeometry *>( fgGeom );
 588   return rv;
 589 }
 590
 591 //______________________________________________________________________
 592 AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
 593                                                 const Text_t* title){
 594     // Returns the pointer of the unique instance
 595
 596     AliEMCALGeometry * rv = 0;
 597     if ( fgGeom == 0 ) {
 598       if ( strcmp(name,"") == 0 ) { // get default geometry
 599          fgGeom = new AliEMCALGeometry(fgDefaultGeometryName, title);
 600       } else {
 601          fgGeom = new AliEMCALGeometry(name, title);
 602       }  // end if strcmp(name,"")
 603       if ( fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
 604       else {
 605          rv = 0;
 606          delete fgGeom;
 607          fgGeom = 0;
 608       } // end if fgInit
 609     }else{
 610         if ( strcmp(fgGeom->GetName(), name) != 0) {
 611           printf("\ncurrent geometry is %s : ", fgGeom->GetName());
 612           printf(" you cannot call %s ",name);
 613         }else{
 614           rv = (AliEMCALGeometry *) fgGeom;
 615         } // end
 616     }  // end if fgGeom
 617     return rv;
 618 }
 619
 620 //_____________________________________________________________________________
 621 Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const {
 622   // Checks whether point is inside the EMCal volume, used in AliEMCALv*.cxx
 623   //
 624   // Code uses cylindrical approximation made of inner radius (for speed)
 625   //
 626   // Points behind EMCAl, i.e. R > outer radius, but eta, phi in acceptance
 627   // are considered to inside
 628
 629   Double_t r=sqrt(x*x+y*y);
 630
 631   if ( r > fEnvelop[0] ) {
 632      Double_t theta;
 633      theta  =    TMath::ATan2(r,z);
 634      Double_t eta;
 635      if(theta == 0)
 636        eta = 9999;
 637      else
 638        eta    =   -TMath::Log(TMath::Tan(theta/2.));
 639      if (eta < fArm1EtaMin || eta > fArm1EtaMax)
 640        return 0;
 641
 642      Double_t phi = TMath::ATan2(y,x) * 180./TMath::Pi();
 643      if (phi < 0) phi += 360;  // phi should go from 0 to 360 in this case
 644      if (phi > fArm1PhiMin && phi < fArm1PhiMax)
 645        return 1;
 646   }
 647   return 0;
 648 }
 649
 650 //
 651 // == Shish-kebab cases ==
 652 //
 653 //________________________________________________________________________________________________
 654 Int_t AliEMCALGeometry::GetAbsCellId(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta) const
 655 {
 656   // 27-aug-04;
 657   // corr. 21-sep-04;
 658   //       13-oct-05; 110 degree case
 659   // May 31, 2006; ALICE numbering scheme:
 660   // 0 <= nSupMod < fNumberOfSuperModules
 661   // 0 <= nModule  < fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1)
 662   // 0 <= nIphi   < fNPHIdiv
 663   // 0 <= nIeta   < fNETAdiv
 664   // 0 <= absid   < fNCells
 665   static Int_t id=0; // have to change from 0 to fNCells-1
 666   if(fKey110DEG == 1 && nSupMod >= 10) { // 110 degree case; last two supermodules
 667     id  = fNCellsInSupMod*10 + (fNCellsInSupMod/2)*(nSupMod-10);
 668   } else {
 669     id  = fNCellsInSupMod*nSupMod;
 670   }
 671   id += fNCellsInModule *nModule;
 672   id += fNPHIdiv *nIphi;
 673   id += nIeta;
 674   if(id<0 || id >= fNCells) {
 675 //     printf(" wrong numerations !!\n");
 676 //     printf("    id      %6i(will be force to -1)\n", id);
 677 //     printf("    fNCells %6i\n", fNCells);
 678 //     printf("    nSupMod %6i\n", nSupMod);
 679 //     printf("    nModule  %6i\n", nModule);
 680 //     printf("    nIphi   %6i\n", nIphi);
 681 //     printf("    nIeta   %6i\n", nIeta);
 682     id = -TMath::Abs(id); // if negative something wrong
 683   }
 684   return id;
 685 }
 686
 687 //________________________________________________________________________________________________
 688 Bool_t  AliEMCALGeometry::CheckAbsCellId(Int_t absId) const
 689 {
 690   // May 31, 2006; only trd1 now
 691   if(absId<0 || absId >= fNCells) return kFALSE;
 692   else                            return kTRUE;
 693 }
 694
 695 //________________________________________________________________________________________________
 696 Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nModule,Int_t &nIphi,Int_t &nIeta) const
 697 {
 698   // 21-sep-04; 19-oct-05;
 699   // May 31, 2006; ALICE numbering scheme:
 700   //
 701   // In:
 702   // absId   - cell is as in Geant,     0<= absId   < fNCells;
 703   // Out:
 704   // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
 705   // nModule  - module number in SM,     0<= nModule  < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
 706   // nIphi   - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv;
 707   // nIeta   - cell number in eta driection inside module; 0<= nIeta < fNETAdiv;
 708   //
 709   static Int_t tmp=0, sm10=0;
 710   if(!CheckAbsCellId(absId)) return kFALSE;
 711
 712   sm10 = fNCellsInSupMod*10;
 713   if(fKey110DEG == 1 && absId >= sm10) { // 110 degree case; last two supermodules
 714     nSupMod = (absId-sm10) / (fNCellsInSupMod/2) + 10;
 715     tmp     = (absId-sm10) % (fNCellsInSupMod/2);
 716   } else {
 717     nSupMod = absId / fNCellsInSupMod;
 718     tmp     = absId % fNCellsInSupMod;
 719   }
 720
 721   nModule  = tmp / fNCellsInModule;
 722   tmp     = tmp % fNCellsInModule;
 723   nIphi   = tmp / fNPHIdiv;
 724   nIeta   = tmp % fNPHIdiv;
 725
 726   return kTRUE;
 727 }
 728
 729 //________________________________________________________________________________________________
 730 void AliEMCALGeometry::GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule,  int &iphim, int &ietam) const
 731 {
 732   // added nSupMod; - 19-oct-05 !
 733   // Alice numbering scheme        - Jun 01,2006
 734   // ietam, iphi - indexes of module in two dimensional grid of SM
 735   // ietam - have to change from 0 to fNZ-1
 736   // iphim - have to change from 0 to nphi-1 (fNPhi-1 or fNPhi/2-1)
 737   static Int_t nphi;
 738
 739   if(fKey110DEG == 1 && nSupMod>=10) nphi = fNPhi/2;
 740   else                               nphi = fNPhi;
 741
 742   ietam = nModule/nphi;
 743   iphim = nModule%nphi;
 744 }
 745
 746 //________________________________________________________________________________________________
 747 void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta,
 748 int &iphi, int &ieta) const
 749 {
 750   //
 751   // Added nSupMod; Nov 25, 05
 752   // Alice numbering scheme  - Jun 01,2006
 753   // IN:
 754   // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
 755   // nModule  - module number in SM,     0<= nModule  < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
 756   // nIphi   - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv;
 757   // nIeta   - cell number in eta driection inside module; 0<= nIeta < fNETAdiv;
 758   //
 759  // OUT:
 760   // ieta, iphi - indexes of cell(tower) in two dimensional grid of SM
 761   // ieta - have to change from 0 to (fNZ*fNETAdiv-1)
 762   // iphi - have to change from 0 to (fNPhi*fNPHIdiv-1 or fNPhi*fNPHIdiv/2-1)
 763   //
 764   static Int_t iphim, ietam;
 765
 766   GetModulePhiEtaIndexInSModule(nSupMod,nModule, iphim, ietam);
 767   //  ieta  = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM)
 768   ieta  = ietam*fNETAdiv + (fNETAdiv - 1 - nIeta); // x(module) = -z(SM)
 769   iphi  = iphim*fNPHIdiv + nIphi;     // y(module) =  y(SM)
 770
 771   if(iphi<0 || ieta<0)
 772   AliDebug(1,Form(" nSupMod %i nModule %i nIphi %i nIeta %i => ieta %i iphi %i\n",
 773   nSupMod, nModule, nIphi, nIeta, ieta, iphi));
 774 }
 775
 776 //________________________________________________________________________________________________
 777 Int_t  AliEMCALGeometry::GetSuperModuleNumber(Int_t absId)  const
 778 {
 779   // Return the number of the  supermodule given the absolute
 780   // ALICE numbering id
 781
 782   static Int_t nSupMod, nModule, nIphi, nIeta;
 783   GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
 784   return nSupMod;
 785 }
 786
 787 //________________________________________________________________________________________________
 788 void  AliEMCALGeometry::GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta,
 789                         Int_t &iphim, Int_t &ietam, Int_t &nModule) const
 790 {
 791   // Transition from cell indexes (ieta,iphi) to module indexes (ietam,iphim, nModule)
 792   static Int_t nphi;
 793   nphi  = GetNumberOfModuleInPhiDirection(nSupMod);
 794
 795   ietam  = ieta/fNETAdiv;
 796   iphim  = iphi/fNPHIdiv;
 797   nModule = ietam * nphi + iphim;
 798 }
 799
 800 //________________________________________________________________________________________________
 801 Int_t  AliEMCALGeometry::GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const
 802 {
 803   // Transition from super module number(nSupMod) and cell indexes (ieta,iphi) to absId
 804   static Int_t ietam, iphim, nModule;
 805   static Int_t nIeta, nIphi; // cell indexes in module
 806
 807   GetModuleIndexesFromCellIndexesInSModule(nSupMod, iphi, ieta, ietam, iphim, nModule);
 808
 809   nIeta = ieta%fNETAdiv;
 810   nIeta = fNETAdiv - 1 - nIeta;
 811   nIphi = iphi%fNPHIdiv;
 812
 813   return GetAbsCellId(nSupMod, nModule, nIphi, nIeta);
 814 }
 815
 816
 817 // Methods for AliEMCALRecPoint - Feb 19, 2006
 818 //________________________________________________________________________________________________
 819 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const
 820 {
 821   // Look to see what the relative
 822   // position inside a given cell is
 823   // for a recpoint.
 824   // Alice numbering scheme - Jun 08, 2006
 825   // In:
 826   // absId   - cell is as in Geant,     0<= absId   < fNCells;
 827   // OUT:
 828   // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
 829
 830   // Shift index taking into account the difference between standard SM
 831   // and SM of half size in phi direction
 832   const Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
 833   static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
 834   if(!CheckAbsCellId(absId)) return kFALSE;
 835
 836   GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
 837   GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
 838
 839   xr = fCentersOfCellsXDir.At(ieta);
 840   zr = fCentersOfCellsEtaDir.At(ieta);
 841
 842   if(nSupMod<10) {
 843     yr = fCentersOfCellsPhiDir.At(iphi);
 844   } else {
 845     yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
 846   }
 847   AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
 848
 849   return kTRUE;
 850 }
 851
 852 //________________________________________________________________________________________________
 853 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t loc[3]) const
 854 {
 855   // Alice numbering scheme - Jun 03, 2006
 856   loc[0] = loc[1] = loc[2]=0.0;
 857   if(RelPosCellInSModule(absId, loc[0],loc[1],loc[2])) {
 858     return kTRUE;
 859   }
 860   return kFALSE;
 861 }
 862
 863 //________________________________________________________________________________________________
 864 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, TVector3 &vloc) const
 865 {
 866   static Double_t loc[3];
 867   if(RelPosCellInSModule(absId,loc)) {
 868     vloc.SetXYZ(loc[0], loc[1], loc[2]);
 869     return kTRUE;
 870   } else {
 871     vloc.SetXYZ(0,0,0);
 872     return kFALSE;
 873   }
 874   // Alice numbering scheme - Jun 03, 2006
 875 }
 876
 877 //________________________________________________________________________________________________
 878 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t distEff, Double_t &xr, Double_t &yr, Double_t &zr) const
 879 {
 880   // Jul 30, 2007 - taking into account position of shower max
 881   // Look to see what the relative
 882   // position inside a given cell is
 883   // for a recpoint.
 884   // In:
 885   // absId   - cell is as in Geant,     0<= absId   < fNCells;
 886   // e       - cluster energy
 887   // OUT:
 888   // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
 889
 890   // Shift index taking into account the difference between standard SM
 891   // and SM of half size in phi direction
 892   const  Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
 893   static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
 894   static Int_t iphim, ietam;
 895   static AliEMCALShishKebabTrd1Module *mod = 0;
 896   static TVector2 v;
 897   if(!CheckAbsCellId(absId)) return kFALSE;
 898
 899   GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
 900   GetModulePhiEtaIndexInSModule(nSupMod, nModule, iphim, ietam);
 901   GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
 902
 903   mod = GetShishKebabModule(ietam);
 904   mod->GetPositionAtCenterCellLine(nIeta, distEff, v);
 905   xr = v.Y() - fParSM[0];
 906   zr = v.X() - fParSM[2];
 907
 908   if(nSupMod<10) {
 909     yr = fCentersOfCellsPhiDir.At(iphi);
 910   } else {
 911     yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
 912   }
 913   AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
 914
 915   return kTRUE;
 916 }
 917
 918 //________________________________________________________________________________________________
 919 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Int_t maxAbsId, Double_t distEff, Double_t &xr, Double_t &yr, Double_t &zr) const
 920 {
 921   // Jul 31, 2007 - taking into account position of shower max and apply coor2.
 922   // Look to see what the relative
 923   // position inside a given cell is
 924   // for a recpoint.
 925   // In:
 926   // absId     - cell is as in Geant,     0<= absId   < fNCells;
 927   // maxAbsId  - abs id of cell with highest energy
 928   // e         - cluster energy
 929   // OUT:
 930   // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
 931
 932   // Shift index taking into account the difference between standard SM
 933   // and SM of half size in phi direction
 934   const  Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
 935   static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
 936   static Int_t iphim, ietam;
 937   static AliEMCALShishKebabTrd1Module *mod = 0;
 938   static TVector2 v;
 939
 940   static Int_t nSupModM, nModuleM, nIphiM, nIetaM, iphiM, ietaM;
 941   static Int_t iphimM, ietamM, maxAbsIdCopy=-1;
 942   static AliEMCALShishKebabTrd1Module *modM = 0;
 943   static Double_t distCorr;
 944
 945   if(!CheckAbsCellId(absId)) return kFALSE;
 946
 947   GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
 948   GetModulePhiEtaIndexInSModule(nSupMod, nModule, iphim, ietam);
 949   GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
 950   mod = GetShishKebabModule(ietam);
 951
 952   if(absId != maxAbsId) {
 953     distCorr = 0.;
 954     if(maxAbsIdCopy != maxAbsId) {
 955       GetCellIndex(maxAbsId, nSupModM, nModuleM, nIphiM, nIetaM);
 956       GetModulePhiEtaIndexInSModule(nSupModM, nModuleM, iphimM, ietamM);
 957       GetCellPhiEtaIndexInSModule(nSupModM,nModuleM,nIphiM,nIetaM, iphiM, ietaM);
 958       modM = GetShishKebabModule(ietamM); // do I need this ?
 959       maxAbsIdCopy = maxAbsId;
 960     }
 961     if(ietamM !=0) {
 962       distCorr = GetEtaModuleSize()*(ietam-ietamM)/TMath::Tan(modM->GetTheta()); // Stay here
 963       //printf(" distCorr %f | dist %f | ietam %i -> etamM %i\n", distCorr, dist, ietam, ietamM);
 964     }
 965     // distEff += distCorr;
 966   }
 967   // Bad resolution in this case, strong bias vs phi
 968   // distEff = 0.0;
 969   mod->GetPositionAtCenterCellLine(nIeta, distEff, v); // Stay here
 970   xr = v.Y() - fParSM[0];
 971   zr = v.X() - fParSM[2];
 972
 973   if(nSupMod<10) {
 974     yr = fCentersOfCellsPhiDir.At(iphi);
 975   } else {
 976     yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
 977   }
 978   AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
 979
 980   return kTRUE;
 981 }
 982
 983 //________________________________________________________________________________________________
 984 void AliEMCALGeometry::CreateListOfTrd1Modules()
 985 {
 986   // Generate the list of Trd1 modules
 987   // which will make up the EMCAL
 988   // geometry
 989
 990   AliDebug(2,Form(" AliEMCALGeometry::CreateListOfTrd1Modules() started "));
 991
 992   AliEMCALShishKebabTrd1Module *mod=0, *mTmp=0; // current module
 993   if(fShishKebabTrd1Modules == 0) {
 994     fShishKebabTrd1Modules = new TList;
 995     fShishKebabTrd1Modules->SetName("ListOfTRD1");
 996     for(int iz=0; iz< GetNZ(); iz++) {
 997       if(iz==0) {
 998         mod  = new AliEMCALShishKebabTrd1Module(TMath::Pi()/2.,this);
 999       } else {
1000         mTmp  = new AliEMCALShishKebabTrd1Module(*mod);
1001         mod   = mTmp;
1002       }
1003       fShishKebabTrd1Modules->Add(mod);
1004     }
1005   } else {
1006     AliDebug(2,Form(" Already exits : "));
1007   }
1008   mod = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(fShishKebabTrd1Modules->GetSize()-1);
1009   fEtaMaxOfTRD1 = mod->GetMaxEtaOfModule(0);
1010
1011   AliDebug(2,Form(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n",
1012                   fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1));
1013   // Feb 20,2006;
1014   // Jun 01, 2006 - ALICE numbering scheme
1015   // define grid for cells in eta(z) and x directions in local coordinates system of SM
1016   // Works just for 2x2 case only -- ?? start here
1017   //
1018   //
1019   // Define grid for cells in phi(y) direction in local coordinates system of SM
1020   // as for 2X2 as for 3X3 - Nov 8,2006
1021   //
1022   AliDebug(2,Form(" Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize()));
1023   Int_t ind=0; // this is phi index
1024   Int_t ieta=0, nModule=0, iphiTemp;
1025   Double_t xr, zr, theta, phi, eta, r, x,y;
1026   TVector3 vglob;
1027   Double_t ytCenterModule=0.0, ytCenterCell=0.0;
1028
1029   fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
1030   fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
1031
1032   Double_t r0 = GetIPDistance() + GetLongModuleSize()/2.;
1033   for(Int_t it=0; it<fNPhi; it++) { // cycle on modules
1034     ytCenterModule = -fParSM[1] + fPhiModuleSize*(2*it+1)/2;  // center of module
1035     for(Int_t ic=0; ic<fNPHIdiv; ic++) { // cycle on cells in module
1036       if(fNPHIdiv==2) {
1037         ytCenterCell = ytCenterModule + fPhiTileSize *(2*ic-1)/2.;
1038       } else if(fNPHIdiv==3){
1039         ytCenterCell = ytCenterModule + fPhiTileSize *(ic-1);
1040       } else if(fNPHIdiv==1){
1041         ytCenterCell = ytCenterModule;
1042       }
1043       fCentersOfCellsPhiDir.AddAt(ytCenterCell,ind);
1044       // Define grid on phi direction
1045       // Grid is not the same for different eta bin;
1046       // Effect is small but is still here
1047       phi = TMath::ATan2(ytCenterCell, r0);
1048       fPhiCentersOfCells.AddAt(phi, ind);
1049
1050       AliDebug(2,Form(" ind %2.2i : y %8.3f ", ind, fCentersOfCellsPhiDir.At(ind)));
1051       ind++;
1052     }
1053   }
1054
1055   fCentersOfCellsEtaDir.Set(fNZ *fNETAdiv);
1056   fCentersOfCellsXDir.Set(fNZ *fNETAdiv);
1057   fEtaCentersOfCells.Set(fNZ *fNETAdiv * fNPhi*fNPHIdiv);
1058   AliDebug(2,Form(" Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize()));
1059   for(Int_t it=0; it<fNZ; it++) {
1060     AliEMCALShishKebabTrd1Module *trd1 = GetShishKebabModule(it);
1061     nModule = fNPhi*it;
1062     for(Int_t ic=0; ic<fNETAdiv; ic++) {
1063       if(fNPHIdiv==2) {
1064         trd1->GetCenterOfCellInLocalCoordinateofSM(ic, xr, zr);      // case of 2X2
1065         GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1066       } if(fNPHIdiv==3) {
1067         trd1->GetCenterOfCellInLocalCoordinateofSM_3X3(ic, xr, zr);  // case of 3X3
1068         GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1069       } if(fNPHIdiv==1) {
1070         trd1->GetCenterOfCellInLocalCoordinateofSM_1X1(xr, zr);      // case of 1X1
1071         GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1072       }
1073       fCentersOfCellsXDir.AddAt(float(xr) - fParSM[0],ieta);
1074       fCentersOfCellsEtaDir.AddAt(float(zr) - fParSM[2],ieta);
1075       // Define grid on eta direction for each bin in phi
1076       for(int iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
1077         x = xr + trd1->GetRadius();
1078         y = fCentersOfCellsPhiDir[iphi];
1079         r = TMath::Sqrt(x*x + y*y + zr*zr);
1080         theta = TMath::ACos(zr/r);
1081         eta   = AliEMCALShishKebabTrd1Module::ThetaToEta(theta);
1082         //        ind   = ieta*fCentersOfCellsPhiDir.GetSize() + iphi;
1083         ind   = iphi*fCentersOfCellsEtaDir.GetSize() + ieta;
1084         fEtaCentersOfCells.AddAt(eta, ind);
1085       }
1086       //printf(" ieta %i : xr + trd1->GetRadius() %f : zr %f : eta %f \n", ieta, xr + trd1->GetRadius(), zr, eta);
1087     }
1088   }
1089   for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1090     AliDebug(2,Form(" ind %2.2i : z %8.3f : x %8.3f", i+1,
1091                     fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i)));
1092   }
1093
1094 }
1095
1096 //________________________________________________________________________________________________
1097 void AliEMCALGeometry::GetGlobal(const Double_t *loc, Double_t *glob, int ind) const
1098 {
1099   // Figure out the global numbering
1100   // of a given supermodule from the
1101   // local numbering and the transformation
1102   // matrix stored by the geometry manager (allows for misaligned
1103   // geometry)
1104
1105   if(ind>=0 && ind < GetNumberOfSuperModules()) {
1106     TString volpath = "ALIC_1/XEN1_1/SMOD_";
1107     volpath += ind+1;
1108
1109     if(GetKey110DEG() && ind>=10) {
1110       volpath = "ALIC_1/XEN1_1/SM10_";
1111       volpath += ind-10+1;
1112     }
1113
1114     if(!gGeoManager->cd(volpath.Data()))
1115       AliFatal(Form("AliEMCALGeometry::GeoManager cannot find path %s!",volpath.Data()));
1116
1117     TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1118     if(m) {
1119       m->LocalToMaster(loc, glob);
1120     } else {
1121       AliFatal("Geo matrixes are not loaded \n") ;
1122     }
1123   }
1124 }
1125
1126 //________________________________________________________________________________________________
1127 void AliEMCALGeometry::GetGlobal(const TVector3 &vloc, TVector3 &vglob, int ind) const
1128 {
1129   //Figure out the global numbering
1130   //of a given supermodule from the
1131   //local numbering given a 3-vector location
1132
1133   static Double_t tglob[3], tloc[3];
1134   vloc.GetXYZ(tloc);
1135   GetGlobal(tloc, tglob, ind);
1136   vglob.SetXYZ(tglob[0], tglob[1], tglob[2]);
1137 }
1138
1139 //________________________________________________________________________________________________
1140 void AliEMCALGeometry::GetGlobal(Int_t absId , double glob[3]) const
1141 {
1142   // Alice numbering scheme - Jun 03, 2006
1143   static Int_t nSupMod, nModule, nIphi, nIeta;
1144   static double loc[3];
1145
1146   if (!gGeoManager || !gGeoManager->IsClosed()) {
1147     AliError("Can't get the global coordinates! gGeoManager doesn't exist or it is still open!");
1148     return;
1149   }
1150
1151   glob[0]=glob[1]=glob[2]=0.0; // bad case
1152   if(RelPosCellInSModule(absId, loc)) {
1153     GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1154
1155     TString volpath = "ALIC_1/XEN1_1/SMOD_";
1156     volpath += (nSupMod+1);
1157
1158     if(GetKey110DEG() && nSupMod>=10) {
1159       volpath = "ALIC_1/XEN1_1/SM10_";
1160       volpath += (nSupMod-10+1);
1161     }
1162     if(!gGeoManager->cd(volpath.Data()))
1163       AliFatal(Form("GeoManager cannot find path %s!",volpath.Data()));
1164
1165     TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1166     if(m) {
1167       m->LocalToMaster(loc, glob);
1168     } else {
1169       AliFatal("Geo matrixes are not loaded \n") ;
1170     }
1171   }
1172 }
1173
1174 //___________________________________________________________________
1175 void AliEMCALGeometry::GetGlobal(Int_t absId , TVector3 &vglob) const
1176 {
1177   // Alice numbering scheme - Jun 03, 2006
1178   static Double_t glob[3];
1179
1180   GetGlobal(absId, glob);
1181   vglob.SetXYZ(glob[0], glob[1], glob[2]);
1182
1183 }
1184
1185 //____________________________________________________________________________
1186 void AliEMCALGeometry::GetGlobal(const AliRecPoint* /*rp*/, TVector3& /* vglob */) const
1187 {
1188   AliFatal(Form("Please use GetGlobalEMCAL(recPoint,gpos) instead of GetGlobal!"));
1189 }
1190
1191 //_________________________________________________________________________________
1192 void AliEMCALGeometry::GetGlobalEMCAL(const AliEMCALRecPoint *rp, TVector3 &vglob) const
1193 {
1194   // Figure out the global numbering
1195   // of a given supermodule from the
1196   // local numbering for RecPoints
1197
1198   static TVector3 vloc;
1199   static Int_t nSupMod, nModule, nIphi, nIeta;
1200
1201   const AliEMCALRecPoint *rpTmp = rp;
1202   const AliEMCALRecPoint *rpEmc = rpTmp;
1203
1204   GetCellIndex(rpEmc->GetAbsId(0), nSupMod, nModule, nIphi, nIeta);
1205   rpTmp->GetLocalPosition(vloc);
1206   GetGlobal(vloc, vglob, nSupMod);
1207 }
1208
1209 //________________________________________________________________________________________________
1210 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Double_t &eta,Double_t &phi) const
1211 {
1212   // Nov 16, 2006- float to double
1213   // version for TRD1 only
1214   static TVector3 vglob;
1215   GetGlobal(absId, vglob);
1216   eta = vglob.Eta();
1217   phi = vglob.Phi();
1218 }
1219
1220 //________________________________________________________________________________________________
1221 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Float_t &eta,Float_t &phi) const
1222 {
1223   // Nov 16,2006 - should be discard in future
1224   static TVector3 vglob;
1225   GetGlobal(absId, vglob);
1226   eta = float(vglob.Eta());
1227   phi = float(vglob.Phi());
1228 }
1229
1230 //________________________________________________________________________________________________
1231 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSM(Int_t nSupMod, Double_t &phiMin, Double_t &phiMax) const
1232 {
1233   // 0<= nSupMod <=11; phi in rad
1234   static int i;
1235   if(nSupMod<0 || nSupMod >11) return kFALSE;
1236   i = nSupMod/2;
1237   phiMin = fPhiBoundariesOfSM[2*i];
1238   phiMax = fPhiBoundariesOfSM[2*i+1];
1239   return kTRUE;
1240 }
1241
1242 //________________________________________________________________________________________________
1243 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const
1244 {
1245   // 0<= nPhiSec <=4; phi in rad
1246   // 0;  gap boundaries between  0th&2th  | 1th&3th SM
1247   // 1;  gap boundaries between  2th&4th  | 3th&5th SM
1248   // 2;  gap boundaries between  4th&6th  | 5th&7th SM
1249   // 3;  gap boundaries between  6th&8th  | 7th&9th SM
1250   // 4;  gap boundaries between  8th&10th | 9th&11th SM
1251   if(nPhiSec<0 || nPhiSec >4) return kFALSE;
1252   phiMin = fPhiBoundariesOfSM[2*nPhiSec+1];
1253   phiMax = fPhiBoundariesOfSM[2*nPhiSec+2];
1254   return kTRUE;
1255 }
1256
1257 //________________________________________________________________________________________________
1258 Bool_t AliEMCALGeometry::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const
1259 {
1260   // Return false if phi belongs a phi cracks between SM
1261
1262   static Int_t i;
1263
1264   if(TMath::Abs(eta) > fEtaMaxOfTRD1) return kFALSE;
1265
1266   phi = TVector2::Phi_0_2pi(phi); // move phi to (0,2pi) boundaries
1267   for(i=0; i<6; i++) {
1268     if(phi>=fPhiBoundariesOfSM[2*i] && phi<=fPhiBoundariesOfSM[2*i+1]) {
1269       nSupMod = 2*i;
1270       if(eta < 0.0) nSupMod++;
1271       AliDebug(1,Form("eta %f phi %f(%5.2f) : nSupMod %i : #bound %i", eta,phi,phi*TMath::RadToDeg(), nSupMod,i));
1272       return kTRUE;
1273     }
1274   }
1275   return kFALSE;
1276 }
1277
1278 //________________________________________________________________________________________________
1279 Bool_t AliEMCALGeometry::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi, Int_t &absId) const
1280 {
1281   // Nov 17,2006
1282   // stay here - phi problem as usual
1283   static Int_t nSupMod, i, ieta, iphi, etaShift, nphi;
1284   static Double_t absEta=0.0, d=0.0, dmin=0.0, phiLoc;
1285   absId = nSupMod = - 1;
1286   if(SuperModuleNumberFromEtaPhi(eta, phi, nSupMod)) {
1287     // phi index first
1288     phi    = TVector2::Phi_0_2pi(phi);
1289     phiLoc = phi - fPhiCentersOfSM[nSupMod/2];
1290     nphi   = fPhiCentersOfCells.GetSize();
1291     if(nSupMod>=10) {
1292       phiLoc = phi - 190.*TMath::DegToRad();
1293       nphi  /= 2;
1294     }
1295
1296     dmin   = TMath::Abs(fPhiCentersOfCells[0]-phiLoc);
1297     iphi   = 0;
1298     for(i=1; i<nphi; i++) {
1299       d = TMath::Abs(fPhiCentersOfCells[i] - phiLoc);
1300       if(d < dmin) {
1301         dmin = d;
1302         iphi = i;
1303       }
1304       //      printf(" i %i : d %f : dmin %f : fPhiCentersOfCells[i] %f \n", i, d, dmin, fPhiCentersOfCells[i]);
1305     }
1306     // odd SM are turned with respect of even SM - reverse indexes
1307     AliDebug(2,Form(" iphi %i : dmin %f (phi %f, phiLoc %f ) ", iphi, dmin, phi, phiLoc));
1308     // eta index
1309     absEta   = TMath::Abs(eta);
1310     etaShift = iphi*fCentersOfCellsEtaDir.GetSize();
1311     dmin     = TMath::Abs(fEtaCentersOfCells[etaShift]-absEta);
1312     ieta     = 0;
1313     for(i=1; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1314       d = TMath::Abs(fEtaCentersOfCells[i+etaShift] - absEta);
1315       if(d < dmin) {
1316         dmin = d;
1317         ieta = i;
1318       }
1319     }
1320     AliDebug(2,Form(" ieta %i : dmin %f (eta=%f) : nSupMod %i ", ieta, dmin, eta, nSupMod));
1321
1322     if(eta<0) iphi = (nphi-1) - iphi;
1323     absId = GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);
1324
1325     return kTRUE;
1326   }
1327   return kFALSE;
1328 }
1329
1330 //________________________________________________________________________________________________
1331 AliEMCALShishKebabTrd1Module* AliEMCALGeometry::GetShishKebabModule(Int_t neta) const
1332 {
1333   //This method was too long to be
1334   //included in the header file - the
1335   //rule checker complained about it's
1336   //length, so we move it here.  It returns the
1337   //shishkebabmodule at a given eta index point.
1338
1339   static AliEMCALShishKebabTrd1Module* trd1=0;
1340   if(fShishKebabTrd1Modules && neta>=0 && neta<fShishKebabTrd1Modules->GetSize()) {
1341     trd1 = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(neta);
1342   } else trd1 = 0;
1343   return trd1;
1344 }
1345
1346 //________________________________________________________________________________________________
1347 Int_t AliEMCALGeometry::GetAbsTRUNumberFromNumberInSm(const Int_t row, const Int_t col, const Int_t sm)
1348 { // Nov 6, 2007
1349   Int_t itru = row + col*GetNModulesInTRUPhi() + sm*GetNTRU();
1350   // printf("  GetAbsTRUNumberFromNumberInSm : row %2i col %2i sm %2i -> itru %2i\n", row, col, sm, itru);
1351   return itru;
1352 }
1353
1354 //________________________________________________________________________________________________
1355 void AliEMCALGeometry::Browse(TBrowser* b)
1356 {
1357   //Browse the modules
1358   if(fShishKebabTrd1Modules) b->Add(fShishKebabTrd1Modules);
1359 }
1360
1361 //________________________________________________________________________________________________
1362 Bool_t AliEMCALGeometry::IsFolder() const
1363 {
1364   //Check if fShishKebabTrd1Modules is in folder
1365   if(fShishKebabTrd1Modules) return kTRUE;
1366   else                       return kFALSE;
1367 }
1368
1369 //________________________________________________________________________________________________
1370 Double_t AliEMCALGeometry::GetPhiCenterOfSM(Int_t nsupmod) const
1371 {
1372   //returns center of supermodule in phi
1373   int i = nsupmod/2;
1374   return fPhiCentersOfSM[i];
1375
1376 }
1377 //____________________________________________________________________________
1378 Bool_t  AliEMCALGeometry::Impact(const TParticle * particle) const
1379 {
1380   // Tells if a particle enters EMCAL
1381   Bool_t in=kFALSE;
1382   Int_t AbsID=0;
1383   TVector3 vtx(particle->Vx(),particle->Vy(),particle->Vz());
1384   TVector3 vimpact(0,0,0);
1385   ImpactOnEmcal(vtx,particle->Theta(),particle->Phi(),AbsID,vimpact);
1386   if(AbsID!=0)
1387     in=kTRUE;
1388   return in;
1389 }
1390 //____________________________________________________________________________
1391 void AliEMCALGeometry::ImpactOnEmcal(TVector3 vtx, Double_t theta, Double_t phi,
1392                                      Int_t & absId, TVector3 & vimpact) const
1393 {
1394   // calculates the impact coordinates on EMCAL (centre of a tower/not on EMCAL surface)
1395   // of a neutral particle
1396   // emitted in the vertex vtx[3] with direction theta and phi in the ALICE global coordinate system
1397
1398   TVector3 p(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) ;
1399
1400   vimpact.SetXYZ(0,0,0);
1401   absId=-1;
1402   if(phi==0 || theta==0) return;
1403
1404    TVector3 direction;
1405    Double_t factor = (GetIPDistance()-vtx[1])/p[1];
1406   direction = vtx + factor*p;
1407
1408   if (!gGeoManager){
1409     AliFatal("Geo manager not initialized\n");
1410   }
1411   //from particle direction -> tower hitted
1412   GetAbsCellIdFromEtaPhi(direction.Eta(),direction.Phi(),absId);
1413
1414   //tower absID hitted -> tower/module plane (evaluated at the center of the tower)
1415   Int_t nSupMod, nModule, nIphi, nIeta;
1416   Double_t loc[3],loc2[3],loc3[3];
1417   Double_t glob[3]={},glob2[3]={},glob3[3]={};
1418
1419   if(!RelPosCellInSModule(absId,loc)) return;
1420
1421   //loc is cell center of tower
1422   GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1423
1424   //look at 2 neighbours-s cell using nIphi={0,1} and nIeta={0,1}
1425   Int_t nIphi2,nIeta2,absId2,absId3;
1426   if(nIeta==0) nIeta2=1;
1427   else nIeta2=0;
1428   absId2=GetAbsCellId(nSupMod,nModule,nIphi,nIeta2);
1429   if(nIphi==0) nIphi2=1;
1430   else nIphi2=0;
1431   absId3=GetAbsCellId(nSupMod,nModule,nIphi2,nIeta);
1432
1433   //2nd point on emcal cell plane
1434   if(!RelPosCellInSModule(absId2,loc2)) return;
1435
1436   //3rd point on emcal cell plane
1437   if(!RelPosCellInSModule(absId3,loc3)) return;
1438
1439   TString volpath = "ALIC_1/XEN1_1/SMOD_";
1440   volpath += (nSupMod+1);
1441
1442   if(GetKey110DEG() && nSupMod>=10) {
1443     volpath = "ALIC_1/XEN1_1/SM10_";
1444     volpath += (nSupMod-10+1);
1445   }
1446   if(!gGeoManager->cd(volpath.Data())){
1447     AliFatal(Form("GeoManager cannot find path %s!",volpath.Data()))
1448     return;
1449   }
1450   TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1451   if(m) {
1452     m->LocalToMaster(loc, glob);
1453     m->LocalToMaster(loc2, glob2);
1454     m->LocalToMaster(loc3, glob3);
1455   } else {
1456     AliFatal("Geo matrixes are not loaded \n") ;
1457   }
1458
1459   //Equation of Plane from glob,glob2,glob3 (Ax+By+Cz+D=0)
1460    Double_t A = glob[1]*(glob2[2]-glob3[2]) + glob2[1]*(glob3[2]-glob[2]) + glob3[1]*(glob[2]-glob2[2]);
1461    Double_t B = glob[2]*(glob2[0]-glob3[0]) + glob2[2]*(glob3[0]-glob[0]) + glob3[2]*(glob[0]-glob2[0]);
1462    Double_t C = glob[0]*(glob2[1]-glob3[1]) + glob2[0]*(glob3[1]-glob[1]) + glob3[0]*(glob[1]-glob2[1]);
1463    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]);
1464   D=-D;
1465
1466   //shift equation of plane from tower/module center to surface along vector (A,B,C) normal to tower/module plane
1467   Double_t dist = GetLongModuleSize()/2.;
1468   Double_t norm = TMath::Sqrt(A*A+B*B+C*C);
1469   Double_t glob4[3]={};
1470   TVector3 dir(A,B,C);
1471   TVector3 point(glob[0],glob[1],glob[2]);
1472   if(point.Dot(dir)<0) dist*=-1;
1473   glob4[0]=glob[0]-dist*A/norm;
1474   glob4[1]=glob[1]-dist*B/norm;
1475   glob4[2]=glob[2]-dist*C/norm;
1476   D = glob4[0]*A +  glob4[1]*B +  glob4[2]*C ;
1477   D = -D;
1478
1479   //Line determination (2 points for equation of line : vtx and direction)
1480   //impact between line (particle) and plane (module/tower plane)
1481   Double_t den = A*(vtx(0)-direction(0)) + B*(vtx(1)-direction(1)) + C*(vtx(2)-direction(2));
1482   if(den==0){
1483     printf("ImpactOnEmcal() No solution :\n");
1484     return;
1485   }
1486
1487   Double_t length = A*vtx(0)+B*vtx(1)+C*vtx(2)+D;
1488   length /=den;
1489
1490   vimpact.SetXYZ(vtx(0)+length*(direction(0)-vtx(0)),vtx(1)+length*(direction(1)-vtx(1)),vtx(2)+length*(direction(2)-vtx(2)));
1491
1492   //shift vimpact from tower/module surface to center along vector (A,B,C) normal to tower/module plane
1493   vimpact.SetXYZ(vimpact(0)+dist*A/norm,vimpact(1)+dist*B/norm,vimpact(2)+dist*C/norm);
1494
1495   return;
1496 }