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