037ea74d297a831241947d58afd100d454dacf2b
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALGeometry.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$*/
17
18 //_________________________________________________________________________
19 // Geometry class  for EMCAL : singleton  
20 // EMCAL consists of layers of scintillator and lead
21 // Places the the Barrel Geometry of The EMCAL at Midrapidity
22 // between 0 and 120 degrees of Phi and
23 // -0.7 to 0.7 in eta 
24 // Number of Modules and Layers may be controlled by 
25 // the name of the instance defined               
26 //*-- Author: Sahal Yacoob (LBL / UCT)
27 //     and  : Yves Schutz (SUBATECH)
28 //     and  : Jennifer Klay (LBL)
29 //     SHASHLYK : Aleksei Pavlinov (WSU)
30
31 // --- AliRoot header files ---
32 #include <TMath.h>
33 #include <TVector3.h>
34
35 // -- ALICE Headers.
36 //#include "AliConst.h"
37
38 // --- EMCAL headers
39 #include "AliEMCALGeometry.h"
40
41 ClassImp(AliEMCALGeometry)
42
43 AliEMCALGeometry *AliEMCALGeometry::fgGeom = 0;
44 Bool_t            AliEMCALGeometry::fgInit = kFALSE;
45 TString name; // contains name of geometry
46
47 //______________________________________________________________________
48 AliEMCALGeometry::~AliEMCALGeometry(void){
49     // dtor
50 }
51
52 //______________________________________________________________________
53 Bool_t AliEMCALGeometry::AreInSameTower(Int_t id1, Int_t id2) const {
54   // Find out whether two hits are in the same tower
55   Int_t idmax = TMath::Max(id1, id2) ; 
56   Int_t idmin = TMath::Min(id1, id2) ;
57   if ( ((idmax - GetNZ() * GetNPhi()) == idmin ) || 
58        ((idmax - 2 * GetNZ() * GetNPhi()) == idmin ) )
59     return kTRUE ; 
60   else 
61     return kFALSE ; 
62 }
63
64 //______________________________________________________________________
65 void AliEMCALGeometry::Init(void){
66   // Initializes the EMCAL parameters
67   // naming convention : GUV_WX_N_ gives the composition of a tower
68   // WX inform about the composition of the EM calorimeter section: 
69   //   thickness in mm of Pb radiator (W) and of scintillator (X), and number of scintillator layers (N)
70   // New geometry: EMCAL_55_25
71   // 24-aug-04 for shish-kebab
72   // SHISH_25 or SHISH_62
73   fgInit = kFALSE; // Assume failed until proven otherwise.
74   name   = GetName();
75   name.ToUpper(); 
76
77   fNZ             = 114;        // granularity along Z (eta) 
78   fNPhi           = 168;        // granularity in phi (azimuth)
79   fArm1PhiMin     = 60.0;       // degrees, Starting EMCAL Phi position
80   fArm1PhiMax     = 180.0;      // degrees, Ending EMCAL Phi position
81   fArm1EtaMin     = -0.7;       // pseudorapidity, Starting EMCAL Eta position
82   fArm1EtaMax     = +0.7;       // pseudorapidity, Ending EMCAL Eta position
83   fIPDistance     = 454.0;      // cm, Radial distance to inner surface of EMCAL
84   fPhiGapForSM    = 0.;         // cm, only for final TRD1 geometry
85
86   // geometry
87   if (name == "EMCAL_55_25") {
88     fECPbRadThickness  = 0.5;  // cm, Thickness of the Pb radiators
89     fECScintThick      = 0.5;  // cm, Thickness of the scintillator
90     fNECLayers         = 25;   // number of scintillator layers
91     
92     fSampling          = 13.1;  // calculated with Birk's law implementation
93  
94     fAlFrontThick      = 3.5;  // cm, Thickness of front Al layer
95     fGap2Active        = 1.0;  // cm, Gap between Al and 1st Scintillator
96   }
97   else if( name == "G56_2_55_19" || name == "EMCAL_5655_21" || name == "G56_2_55_19_104_14"|| name == "G65_2_64_19" || name == "EMCAL_6564_21"){
98     Fatal("Init", "%s is an old geometry! Please update your Config file", name.Data()) ;
99   }
100   else if(name.Contains("SHISH")){
101     // 7-sep-05; integration issue
102     fArm1PhiMin     = 80.0;     // 60  -> 80
103     fArm1PhiMax     = 180.0;    // 180 -> 200
104
105     fNumberOfSuperModules = 10; // 12 = 6 * 2 (6 in phi, 2 in Z);
106     fSteelFrontThick = 2.54;    //  9-sep-04
107     fIPDistance      = 460.0;
108     fFrontSteelStrip = fPassiveScintThick = 0.0; // 13-may-05
109     fLateralSteelStrip = 0.025; // before MAY 2005 
110     fPhiModuleSize   = fEtaModuleSize   = 11.4;
111     fPhiTileSize = fEtaTileSize      = 5.52; // (11.4-5.52*2)/2. = 0.18 cm (wall thickness)
112     fNPhi            = 14;
113     fNZ              = 30;
114     fAlFrontThick    = fGap2Active = 0;
115     fNPHIdiv = fNETAdiv = 2;
116
117     fNECLayers       = 62;
118     fECScintThick    = fECPbRadThickness = 0.2;
119     fSampling        = 1.;  // 30-aug-04 - should be calculated
120     if(name.Contains("TWIST")) { // all about EMCAL module
121       fNZ             = 27;  // 16-sep-04
122     } else if(name.Contains("TRD")) {
123       fIPDistance      = 428.0;  //  11-may-05
124       fSteelFrontThick = 0.0;    // 3.17 -> 0.0; 28-mar-05 : no stell plate
125       fNPhi            = 12;
126       fSampling       = 12.327;
127       fPhiModuleSize = fEtaModuleSize = 12.26;
128       fNZ            = 26;     // 11-oct-04
129       fTrd1Angle     = 1.3;    // in degree
130 // 18-nov-04; 1./0.08112=12.327
131 // http://pdsfweb01.nersc.gov/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
132       if(name.Contains("TRD1")) {       // 30-jan-05
133         // for final design
134         fPhiGapForSM    = 2.;         // cm, only for final TRD1 geometry
135         if(name.Contains("MAY05") || name.Contains("WSUC") || name.Contains("FINAL")){
136           fNumberOfSuperModules = 12; // 20-may-05
137           if(name.Contains("WSUC")) fNumberOfSuperModules = 1; // 27-may-05
138           fNECLayers     = 77;       // (13-may-05 from V.Petrov)
139           fPhiModuleSize = 12.5;     // 20-may-05 - rectangular shape
140           fEtaModuleSize = 11.9;
141           fECScintThick  = fECPbRadThickness = 0.16;// (13-may-05 from V.Petrov)
142           fFrontSteelStrip   = 0.025;// 0.025cm = 0.25mm  (13-may-05 from V.Petrov)
143           fLateralSteelStrip = 0.01; // 0.01cm  = 0.1mm   (13-may-05 from V.Petrov) - was 0.025
144           fPassiveScintThick = 0.8;  // 0.8cm   = 8mm     (13-may-05 from V.Petrov)
145           fNZ                = 24;
146           fTrd1Angle         = 1.5;  // 1.3 or 1.5
147
148           if(name.Contains("FINAL")) { // 9-sep-05
149             fNumberOfSuperModules = 10;
150             fPhiModuleSize = 12.26 - fPhiGapForSM / Float_t(fNPhi); // first assumption
151             fEtaModuleSize = fPhiModuleSize;
152           }
153         }
154       } else if(name.Contains("TRD2")) {       // 30-jan-05
155         fSteelFrontThick = 0.0;         // 11-mar-05
156         fIPDistance+= fSteelFrontThick; // 1-feb-05 - compensate absence of steel plate
157         fTrd1Angle  = 1.64;             // 1.3->1.64
158         fTrd2AngleY = fTrd1Angle;       //  symmetric case now
159         fEmptySpace    = 0.2; // 2 mm
160         fTubsR         = fIPDistance; // 31-jan-05 - as for Fred case
161
162         fPhiModuleSize  = fTubsR*2.*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.);
163         fPhiModuleSize -= fEmptySpace/2.; // 11-mar-05  
164         fEtaModuleSize  = fPhiModuleSize; // 20-may-05 
165         fTubsTurnAngle  = 3.;
166       }
167       fNPHIdiv = fNETAdiv  = 2;   // 13-oct-04 - division again
168       if(name.Contains("3X3")) {   // 23-nov-04
169         fNPHIdiv = fNETAdiv  = 3;
170       } else if(name.Contains("4X4")) {
171         fNPHIdiv = fNETAdiv  = 4;
172       }
173     }
174     fPhiTileSize = fPhiModuleSize/2. - fLateralSteelStrip; // 13-may-05 
175     fEtaTileSize = fEtaModuleSize/2. - fLateralSteelStrip; // 13-may-05 
176
177     if(name.Contains("25")){
178       fNECLayers     = 25;
179       fECScintThick  = fECPbRadThickness = 0.5;
180     }
181     if(name.Contains("WSUC")){ // 18-may-05 - about common structure
182       fShellThickness = 30.; // should be change 
183       fNPhi = fNZ = 4; 
184     }
185     // constant for transition absid <--> indexes
186     fNCellsInTower  = fNPHIdiv*fNETAdiv;
187     fNCellsInSupMod = fNCellsInTower*fNPhi*fNZ;
188     fNCells         = fNCellsInSupMod*fNumberOfSuperModules;
189
190     fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness);
191     if(name.Contains("MAY05")) fLongModuleSize += (fFrontSteelStrip + fPassiveScintThick);
192
193     // 30-sep-04
194     if(name.Contains("TRD")) {
195       f2Trd1Dx2 = fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.);
196       if(name.Contains("TRD2")) {  // 27-jan-05
197         f2Trd2Dy2 = fPhiModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.);
198       }
199     }
200   }
201   else
202     Fatal("Init", "%s is an undefined geometry!", name.Data()) ; 
203                  
204
205   fNPhiSuperModule = fNumberOfSuperModules/2;
206   if(fNPhiSuperModule<1) fNPhiSuperModule = 1;
207   //There is always one more scintillator than radiator layer because of the first block of aluminium
208   fShellThickness = fAlFrontThick + fGap2Active + fNECLayers*GetECScintThick()+(fNECLayers-1)*GetECPbRadThick();
209   if(name.Contains("SHISH")) {
210     fShellThickness = fSteelFrontThick + fLongModuleSize;
211     if(name.Contains("TWIST")) { // 13-sep-04
212       fShellThickness  = TMath::Sqrt(fLongModuleSize*fLongModuleSize + fPhiModuleSize*fEtaModuleSize);
213       fShellThickness += fSteelFrontThick;
214     } else if(name.Contains("TRD")) { // 1-oct-04
215       fShellThickness  = TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2);
216       fShellThickness += fSteelFrontThick;
217     }
218   }
219
220   fZLength        = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax); // Z coverage
221   fEnvelop[0]     = fIPDistance; // mother volume inner radius
222   fEnvelop[1]     = fIPDistance + fShellThickness; // mother volume outer r.
223   fEnvelop[2]     = 1.00001*fZLength; // add some padding for mother volume. 
224   
225   fgInit = kTRUE; 
226   
227   if (kTRUE) {
228     printf("Init: geometry of EMCAL named %s is as follows:\n", name.Data());
229     printf( "               ECAL      : %d x (%f cm Pb, %f cm Sc) \n", GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ; 
230     if(name.Contains("SHISH")){
231       printf(" fIPDistance       %6.3f cm \n", fIPDistance);
232       if(fSteelFrontThick>0.) 
233       printf(" fSteelFrontThick  %6.3f cm \n", fSteelFrontThick);
234       printf(" fNPhi %i   |  fNZ %i \n", fNPhi, fNZ);
235       if(name.Contains("MAY05")){
236         printf(" fFrontSteelStrip         %6.4f cm (thickness of front steel strip)\n", 
237         fFrontSteelStrip);
238         printf(" fLateralSteelStrip       %6.4f cm (thickness of lateral steel strip)\n", 
239         fLateralSteelStrip);
240         printf(" fPassiveScintThick  %6.4f cm (thickness of front passive Sc tile)\n",
241         fPassiveScintThick);
242       }
243       printf(" X:Y module size   %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize);
244       printf(" X:Y   tile size   %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize);
245       printf(" fLongModuleSize   %6.3f cm \n", fLongModuleSize);
246       printf(" #supermodule in phi direction %i \n", fNPhiSuperModule );
247     }
248     if(name.Contains("TRD")) {
249       printf(" fTrd1Angle %7.4f\n", fTrd1Angle);
250       printf(" f2Trd1Dx2  %7.4f\n",  f2Trd1Dx2);
251       if(name.Contains("TRD2")) {
252         printf(" fTrd2AngleY     %7.4f\n", fTrd2AngleY);
253         printf(" f2Trd2Dy2       %7.4f\n", f2Trd2Dy2);
254         printf(" fTubsR          %7.2f cm\n", fTubsR);
255         printf(" fTubsTurnAngle  %7.4f\n", fTubsTurnAngle);
256         printf(" fEmptySpace     %7.4f cm\n", fEmptySpace);
257       } else if(name.Contains("TRD1") && name.Contains("FINAL")){
258         printf(" fPhiGapForSM  %7.4f cm \n",  fPhiGapForSM);
259       }
260     }
261     printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
262     printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f\n",  
263            GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() );
264   }
265 }
266
267 //______________________________________________________________________
268 AliEMCALGeometry *  AliEMCALGeometry::GetInstance(){ 
269   // Returns the pointer of the unique instance
270   
271   return static_cast<AliEMCALGeometry *>( fgGeom ) ; 
272 }
273
274 //______________________________________________________________________
275 AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
276                                                 const Text_t* title){
277     // Returns the pointer of the unique instance
278
279     AliEMCALGeometry * rv = 0; 
280     if ( fgGeom == 0 ) {
281         if ( strcmp(name,"") == 0 ) rv = 0;
282         else {    
283             fgGeom = new AliEMCALGeometry(name, title);
284             if ( fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
285             else {
286                 rv = 0; 
287                 delete fgGeom; 
288                 fgGeom = 0; 
289             } // end if fgInit
290         } // end if strcmp(name,"")
291     }else{
292         if ( strcmp(fgGeom->GetName(), name) != 0 ) {
293           printf("\ncurrent geometry is ") ;  
294           printf(fgGeom->GetName());
295           printf("\n                      you cannot call     "); 
296           printf(name);  
297         }else{
298           rv = (AliEMCALGeometry *) fgGeom; 
299         } // end if
300     }  // end if fgGeom
301     return rv; 
302 }
303
304 //______________________________________________________________________
305 Int_t AliEMCALGeometry::TowerIndex(Int_t ieta,Int_t iphi) const {
306   // Returns the tower index number from the based on the Z and Phi
307   // index numbers.
308   // Inputs:
309   //   Int_t ieta    // index along z axis [1-fNZ]
310   //   Int_t iphi  // index along phi axis [1-fNPhi]
311   // Outputs:
312   //   none.
313   // Returned
314   //   Int_t index // Tower index number 
315   
316   if ( (ieta <= 0 || ieta>GetNEta()) || 
317        (iphi <= 0 || iphi>GetNPhi())) {
318     Error("TowerIndex", "Unexpected parameters eta = %d phi = %d!", ieta, iphi) ; 
319     return -1;
320   }
321   return ( (iphi - 1)*GetNEta() + ieta ); 
322 }
323
324 //______________________________________________________________________
325 void AliEMCALGeometry::TowerIndexes(Int_t index,Int_t &ieta,Int_t &iphi) const {
326   // Inputs:
327   //   Int_t index // Tower index number [1-fNZ*fNPhi]
328   // Outputs:
329   //   Int_t ieta    // index allong z axis [1-fNZ]
330   //   Int_t iphi  // index allong phi axis [1-fNPhi]
331   // Returned
332   //   none.
333
334   Int_t nindex = 0;
335
336   if ( IsInECA(index) ) { // ECAL index
337     nindex = index ;
338   }
339   else {
340     Error("TowerIndexes", "Unexpected Id number!") ;
341     ieta = -1;
342     iphi = -1;
343     return;
344   }   
345
346   if (nindex%GetNZ()) 
347     iphi = nindex / GetNZ() + 1 ; 
348   else 
349     iphi = nindex / GetNZ() ; 
350   ieta = nindex - (iphi - 1) * GetNZ() ; 
351
352   if (gDebug==2)
353     printf("TowerIndexes: index=%d,%d, ieta=%d, iphi = %d", index, nindex,ieta, iphi) ; 
354   return;
355   
356 }
357
358 //______________________________________________________________________
359 void AliEMCALGeometry::EtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi) const {
360     // given the tower index number it returns the based on the eta and phi
361     // of the tower.
362     // Inputs:
363     //   Int_t index // Tower index number [1-fNZ*fNPhi]
364     // Outputs:
365     //   Float_t eta  // eta of center of tower in pseudorapidity
366     //   Float_t phi  // phi of center of tower in degrees
367     // Returned
368     //   none.
369     Int_t ieta, iphi;
370     Float_t deta, dphi ;
371
372     TowerIndexes(index,ieta,iphi);
373     
374     if (gDebug == 2) 
375       printf("EtaPhiFromIndex: index = %d, ieta = %d, iphi = %d", index, ieta, iphi) ;
376
377     deta = (GetArm1EtaMax()-GetArm1EtaMin())/(static_cast<Float_t>(GetNEta()));
378     eta  = GetArm1EtaMin() + ((static_cast<Float_t>(ieta) - 0.5 ))*deta;
379
380     dphi = (GetArm1PhiMax() - GetArm1PhiMin())/(static_cast<Float_t>(GetNPhi()));  // in degrees.
381     phi  = GetArm1PhiMin() + dphi*(static_cast<Float_t>(iphi) - 0.5);//iphi range [1-fNphi].
382 }
383
384 //______________________________________________________________________
385 Int_t AliEMCALGeometry::TowerIndexFromEtaPhi(Float_t eta,Float_t phi) const {
386     // returns the tower index number based on the eta and phi of the tower.
387     // Inputs:
388     //   Float_t eta  // eta of center of tower in pseudorapidity
389     //   Float_t phi  // phi of center of tower in degrees
390     // Outputs:
391     //   none.
392     // Returned
393     //   Int_t index // Tower index number [1-fNZ*fNPhi]
394
395     Int_t ieta,iphi;
396
397     ieta = static_cast<Int_t> ( 1 + (static_cast<Float_t>(GetNEta()) * (eta - GetArm1EtaMin()) / (GetArm1EtaMax() - GetArm1EtaMin())) ) ;
398
399     if( ieta <= 0 || ieta > GetNEta() ) { 
400       Error("TowerIndexFromEtaPhi", "Unexpected (eta, phi) = (%f, %f) value, outside of EMCAL!", eta, phi) ; 
401       return -1 ; 
402     }
403
404     iphi = static_cast<Int_t> ( 1 + (static_cast<Float_t>(GetNPhi()) * (phi - GetArm1PhiMin()) / (GetArm1PhiMax() - GetArm1PhiMin())) ) ;
405
406     if( iphi <= 0 || iphi > GetNPhi() ) { 
407       Error("TowerIndexFromEtaPhi", "Unexpected (eta, phi) = (%f, %f) value, outside of EMCAL!", eta, phi) ; 
408       return -1 ; 
409     }
410
411     return TowerIndex(ieta,iphi);
412 }
413
414 //______________________________________________________________________
415 Bool_t AliEMCALGeometry::AbsToRelNumbering(Int_t AbsId, Int_t *relid) const {
416     // Converts the absolute numbering into the following array/
417     //  relid[0] = Row number inside EMCAL
418     //  relid[1] = Column number inside EMCAL
419     // Input:
420     //   Int_t AbsId // Tower index number [1-2*fNZ*fNPhi]
421     // Outputs:
422     //   Int_t *relid // array of 2. Described above.
423     Bool_t rv  = kTRUE ;
424     Int_t ieta=0,iphi=0,index=AbsId;
425
426     TowerIndexes(index,ieta,iphi);
427     relid[0] = ieta;
428     relid[1] = iphi;
429
430     return rv;
431 }
432
433 //______________________________________________________________________
434 void AliEMCALGeometry::PosInAlice(const Int_t *relid, Float_t &theta, Float_t &phi) const 
435 {
436   // Converts the relative numbering into the local EMCAL-module (x, z)
437   // coordinates
438   Int_t ieta = relid[0]; // offset along x axis
439   Int_t iphi = relid[1]; // offset along z axis
440   Int_t index;
441   Float_t eta;
442   
443   index = TowerIndex(ieta,iphi);
444   EtaPhiFromIndex(index,eta,phi);
445   //theta = 180.*(2.0*TMath::ATan(TMath::Exp(-eta)))/TMath::Pi();
446   theta = 2.0*TMath::ATan(TMath::Exp(-eta));
447
448   // correct for distance to IP
449   Float_t d = GetIP2ECASection() - GetIPDistance() ;  
450
451   Float_t correction = 1 + d/GetIPDistance() ; 
452   Float_t tantheta = TMath::Tan(theta) * correction ; 
453   theta = TMath::ATan(tantheta) * TMath::RadToDeg() ; 
454   if (theta < 0 ) 
455     theta += 180. ; 
456   
457   return;
458 }
459
460 //______________________________________________________________________
461 void AliEMCALGeometry::PosInAlice(Int_t absid, Float_t &theta, Float_t &phi) const 
462 {
463   // Converts the relative numbering into the local EMCAL-module (x, z)
464   // coordinates
465   Int_t relid[2] ; 
466   AbsToRelNumbering(absid, relid) ;
467   Int_t ieta = relid[0]; // offset along x axis
468   Int_t iphi = relid[1]; // offset along z axis
469   Int_t index;
470   Float_t eta;
471   
472   index = TowerIndex(ieta,iphi);
473   EtaPhiFromIndex(index,eta,phi);
474   theta = 2.0*TMath::ATan(TMath::Exp(-eta)) ;
475   
476   // correct for distance to IP
477   Float_t d = 0. ; 
478   if (IsInECA(absid))
479     d = GetIP2ECASection() - GetIPDistance() ; 
480   else {
481     Error("PosInAlice", "Unexpected id # %d!", absid) ; 
482     return;
483   }
484
485   Float_t correction = 1 + d/GetIPDistance() ; 
486   Float_t tantheta = TMath::Tan(theta) * correction ; 
487   theta = TMath::ATan(tantheta) * TMath::RadToDeg() ; 
488   if (theta < 0 ) 
489     theta += 180. ; 
490   
491   return;
492 }
493
494 //______________________________________________________________________
495 void AliEMCALGeometry::XYZFromIndex(const Int_t *relid,Float_t &x,Float_t &y, Float_t &z) const {
496     // given the tower relative number it returns the X, Y and Z
497     // of the tower.
498     
499     // Outputs:
500     //   Float_t x  // x of center of tower in cm
501     //   Float_t y  // y of center of tower in cm
502     //   Float_t z  // z of centre of tower in cm
503     // Returned
504     //   none.
505     
506     Float_t eta,theta, phi,cylradius=0. ;
507     
508     Int_t ieta = relid[0]; // offset along x axis
509     Int_t iphi = relid[1]; // offset along z axis.
510     Int_t index;
511     
512     index = TowerIndex(ieta,iphi);
513     EtaPhiFromIndex(index,eta,phi);
514     theta = 180.*(2.0*TMath::ATan(TMath::Exp(-eta)))/TMath::Pi();
515     
516     cylradius = GetIP2ECASection() ;  
517
518     Double_t  kDeg2Rad = TMath::DegToRad() ; 
519     x =  cylradius * TMath::Cos(phi * kDeg2Rad ) ;
520     y =  cylradius * TMath::Sin(phi * kDeg2Rad ) ; 
521     z =  cylradius / TMath::Tan(theta * kDeg2Rad ) ; 
522  
523  return;
524
525
526 //______________________________________________________________________
527 void AliEMCALGeometry::XYZFromIndex(Int_t absid,  TVector3 &v) const {
528     // given the tower relative number it returns the X, Y and Z
529     // of the tower.
530     
531     // Outputs:
532     //   Float_t x  // x of center of tower in cm
533     //   Float_t y  // y of center of tower in cm
534     //   Float_t z  // z of centre of tower in cm
535     // Returned
536     //   none.
537     
538     Float_t theta, phi,cylradius=0. ;
539         
540     PosInAlice(absid, theta, phi) ; 
541     
542     if ( IsInECA(absid) ) 
543       cylradius = GetIP2ECASection() ;
544     else {
545       Error("XYZFromIndex", "Unexpected Tower section") ;
546       return;
547     }
548
549     Double_t  kDeg2Rad = TMath::DegToRad() ; 
550     v.SetX(cylradius * TMath::Cos(phi * kDeg2Rad ) );
551     v.SetY(cylradius * TMath::Sin(phi * kDeg2Rad ) ); 
552     v.SetZ(cylradius / TMath::Tan(theta * kDeg2Rad ) ) ; 
553  
554  return;
555
556
557 Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const {
558   // Checks whether point is inside the EMCal volume
559   //
560   // Code uses cylindrical approximation made of inner radius (for speed)
561   //
562   // Points behind EMCAl, i.e. R > outer radius, but eta, phi in acceptance 
563   // are considered to inside
564
565   Double_t r=sqrt(x*x+y*y);
566
567   if ( r > fEnvelop[0] ) {
568      Double_t theta;
569      theta  =    TMath::ATan2(r,z);
570      Double_t eta;
571      if(theta == 0) 
572        eta = 9999;
573      else 
574        eta    =   -TMath::Log(TMath::Tan(theta/2.));
575      if (eta < fArm1EtaMin || eta > fArm1EtaMax)
576        return 0;
577  
578      Double_t phi = TMath::ATan2(y,x) * 180./TMath::Pi();
579      if (phi > fArm1PhiMin && phi < fArm1PhiMax)
580        return 1;
581   }
582   return 0;
583 }
584
585 //
586 // == Shish-kebab cases ==
587 //
588 Int_t AliEMCALGeometry::GetAbsCellId(const int nSupMod, const int nTower, const int nIphi, const int nIeta)
589 { // 27-aug-04; corr. 21-sep-04
590   static Int_t id; // have to change from 1 to fNCells
591   id  = fNCellsInSupMod*(nSupMod-1);
592   id += fNCellsInTower *(nTower-1);
593   id += fNPHIdiv *(nIphi-1);
594   id += nIeta;
595   if(id<=0 || id > fNCells) {
596 //     printf(" wrong numerations !!\n");
597 //     printf("    id      %6i(will be force to -1)\n", id);
598 //     printf("    fNCells %6i\n", fNCells);
599 //     printf("    nSupMod %6i\n", nSupMod);
600 //     printf("    nTower  %6i\n", nTower);
601 //     printf("    nIphi   %6i\n", nIphi);
602 //     printf("    nIeta   %6i\n", nIeta);
603     id = -1;
604   }
605   return id;
606 }
607
608 Bool_t  AliEMCALGeometry::CheckAbsCellId(Int_t ind)
609 { // 17-niv-04 - analog of IsInECA
610    if(name.Contains("TRD")) {
611      if(ind<=0 || ind > fNCells) return kFALSE;
612      else                        return kTRUE;
613    } else return IsInECA(ind);
614 }
615
616 Bool_t AliEMCALGeometry::GetCellIndex(const Int_t absId,Int_t &nSupMod,Int_t &nTower,Int_t &nIphi,Int_t &nIeta)
617 { // 21-sep-04
618   static Int_t tmp=0;
619   if(absId<=0 || absId>fNCells) {
620 //     Info("GetCellIndex"," wrong abs Id %i !! \n", absId); 
621     return kFALSE;
622   }
623   nSupMod = (absId-1) / fNCellsInSupMod + 1;
624   tmp     = (absId-1) % fNCellsInSupMod;
625
626   nTower  = tmp / fNCellsInTower + 1;
627   tmp     = tmp % fNCellsInTower;
628
629   nIphi     = tmp / fNPHIdiv + 1;
630   nIeta     = tmp % fNPHIdiv + 1;
631
632   return kTRUE;
633 }
634
635 void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(const int nTower, const int nIphi, const int nIeta, 
636 int &iphi, int &ieta)
637 { // don't check validity of nTower, nIphi and nIeta index
638   // have to change  - 1-nov-04 ?? 
639   static Int_t iphit, ietat;
640
641   ietat = (nTower-1)/fNPhi;
642   ieta  = ietat*fNETAdiv + nIeta; // change from 1 to fNZ*fNETAdiv
643
644   iphit = (nTower-1)%fNPhi;
645   iphi  = iphit*fNPHIdiv + nIphi;  // change from 1 to fNPhi*fNPHIdiv
646 }