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