New methods in geometry, update geometry for 3x3 case
[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 80 and 180(or 190) 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 //     EMCAL geometry tree:
27 //     EMCAL -> superModule -> module -> cell
28 //     Indexes
29 //     absId -> nSupMod     -> nTower -> (nIphi,nIeta)
30 //
31 //*-- Author: Sahal Yacoob (LBL / UCT)
32 //     and  : Yves Schutz (SUBATECH)
33 //     and  : Jennifer Klay (LBL)
34 //     SHASHLYK : Aleksei Pavlinov (WSU)
35 //
36 // --- AliRoot header files ---
37 #include <assert.h>
38 #include <Riostream.h>
39
40 #include <TMath.h>
41 #include <TVector3.h>
42 #include <TObjArray.h>
43 #include <TGeoManager.h>
44 #include <TGeoNode.h>
45 #include <TGeoMatrix.h>
46 #include <TMatrixD.h>
47 #include <TObjString.h>
48 #include <TClonesArray.h>
49 #include <TBrowser.h>
50
51 // -- ALICE Headers.
52 //#include "AliConst.h"
53 #include "AliLog.h"
54
55 // --- EMCAL headers
56 #include "AliEMCALGeometry.h"
57 #include "AliEMCALShishKebabTrd1Module.h"
58 #include "AliEMCALRecPoint.h"
59 #include "AliEMCALDigit.h"
60 #include "AliEMCALHistoUtilities.h"
61
62 ClassImp(AliEMCALGeometry)
63
64 // these initialisations are needed for a singleton
65 AliEMCALGeometry  *AliEMCALGeometry::fgGeom      = 0;
66 Bool_t             AliEMCALGeometry::fgInit      = kFALSE;
67
68
69 AliEMCALGeometry::AliEMCALGeometry() 
70   : AliGeometry(),
71     fGeoName(0),fArrayOpts(0),fAlFrontThick(0.),fECPbRadThickness(0.),fECScintThick(0.),
72     fNECLayers(0),fArm1PhiMin(0.),fArm1PhiMax(0.),fArm1EtaMin(0.),fArm1EtaMax(0.),fIPDistance(0.),
73     fShellThickness(0.),fZLength(0.),fGap2Active(0.),fNZ(0),fNPhi(0),fSampling(0.),fNumberOfSuperModules(0),
74     fSteelFrontThick(0.),fFrontSteelStrip(0.),fLateralSteelStrip(0.),fPassiveScintThick(0.),fPhiModuleSize(0.),
75     fEtaModuleSize(0.),fPhiTileSize(0.),fEtaTileSize(0.),fLongModuleSize(0.),fNPhiSuperModule(0),fNPHIdiv(0),fNETAdiv(0),
76     fNCells(0),fNCellsInSupMod(0),fNCellsInTower(0),fNTRU(0),fNTRUEta(0),fNTRUPhi(0),fTrd1Angle(0.),f2Trd1Dx2(0.),
77     fPhiGapForSM(0.),fKey110DEG(0),fPhiBoundariesOfSM(0), fPhiCentersOfSM(0),fEtaMaxOfTRD1(0),
78     fTrd2AngleY(0.),f2Trd2Dy2(0.),fEmptySpace(0.),fTubsR(0.),fTubsTurnAngle(0.),fCentersOfCellsEtaDir(0),
79     fCentersOfCellsXDir(0),fCentersOfCellsPhiDir(0),fEtaCentersOfCells(0),fPhiCentersOfCells(0),
80     fShishKebabTrd1Modules(0),fNAdditionalOpts(0) 
81
82   // default ctor only for internal usage (singleton)
83   // must be kept public for root persistency purposes, but should never be called by the outside world    
84   //  CreateListOfTrd1Modules();
85   AliDebug(2, "AliEMCALGeometry : default ctor ");
86 }
87 //______________________________________________________________________
88 AliEMCALGeometry::AliEMCALGeometry(const Text_t* name, const Text_t* title) 
89   : AliGeometry(name, title),
90     fGeoName(0),fArrayOpts(0),fAlFrontThick(0.),fECPbRadThickness(0.),fECScintThick(0.),
91     fNECLayers(0),fArm1PhiMin(0.),fArm1PhiMax(0.),fArm1EtaMin(0.),fArm1EtaMax(0.),fIPDistance(0.),
92     fShellThickness(0.),fZLength(0.),fGap2Active(0.),fNZ(0),fNPhi(0),fSampling(0.),fNumberOfSuperModules(0),
93     fSteelFrontThick(0.),fFrontSteelStrip(0.),fLateralSteelStrip(0.),fPassiveScintThick(0.),fPhiModuleSize(0.),
94     fEtaModuleSize(0.),fPhiTileSize(0.),fEtaTileSize(0.),fLongModuleSize(0.),fNPhiSuperModule(0),fNPHIdiv(0),fNETAdiv(0),
95     fNCells(0),fNCellsInSupMod(0),fNCellsInTower(0),fNTRU(0),fNTRUEta(0),fNTRUPhi(0),fTrd1Angle(0.),f2Trd1Dx2(0.),
96     fPhiGapForSM(0.),fKey110DEG(0),fPhiBoundariesOfSM(0), fPhiCentersOfSM(0), fEtaMaxOfTRD1(0),
97     fTrd2AngleY(0.),f2Trd2Dy2(0.),fEmptySpace(0.),fTubsR(0.),fTubsTurnAngle(0.),fCentersOfCellsEtaDir(0),
98     fCentersOfCellsXDir(0),fCentersOfCellsPhiDir(0),fEtaCentersOfCells(0),fPhiCentersOfCells(0),
99     fShishKebabTrd1Modules(0),fNAdditionalOpts(0)
100 {
101   // ctor only for internal usage (singleton)
102   AliDebug(2, Form("AliEMCALGeometry(%s,%s) ", name,title));
103
104   Init();
105
106   CreateListOfTrd1Modules();
107
108   if (AliDebugLevel()>=2) {
109     PrintGeometry();
110   }
111
112 }
113 //______________________________________________________________________
114 AliEMCALGeometry::AliEMCALGeometry(const AliEMCALGeometry& geom)
115   : AliGeometry(geom),
116     fGeoName(geom.fGeoName),
117     fArrayOpts(geom.fArrayOpts),
118     fAlFrontThick(geom.fAlFrontThick),
119     fECPbRadThickness(geom.fECPbRadThickness),
120     fECScintThick(geom.fECScintThick),
121     fNECLayers(geom.fNECLayers),
122     fArm1PhiMin(geom.fArm1PhiMin),
123     fArm1PhiMax(geom.fArm1PhiMax),
124     fArm1EtaMin(geom.fArm1EtaMin),
125     fArm1EtaMax(geom.fArm1EtaMax),
126     fIPDistance(geom.fIPDistance),
127     fShellThickness(geom.fShellThickness),
128     fZLength(geom.fZLength),
129     fGap2Active(geom.fGap2Active),
130     fNZ(geom.fNZ),
131     fNPhi(geom.fNPhi),
132     fSampling(geom.fSampling),
133     fNumberOfSuperModules(geom.fNumberOfSuperModules),
134     fSteelFrontThick(geom.fSteelFrontThick),
135     fFrontSteelStrip(geom.fFrontSteelStrip),
136     fLateralSteelStrip(geom.fLateralSteelStrip),
137     fPassiveScintThick(geom.fPassiveScintThick),
138     fPhiModuleSize(geom.fPhiModuleSize),
139     fEtaModuleSize(geom.fEtaModuleSize),
140     fPhiTileSize(geom.fPhiTileSize),
141     fEtaTileSize(geom.fEtaTileSize),
142     fLongModuleSize(geom.fLongModuleSize),
143     fNPhiSuperModule(geom.fNPhiSuperModule),
144     fNPHIdiv(geom.fNPHIdiv),
145     fNETAdiv(geom.fNETAdiv),
146     fNCells(geom.fNCells),
147     fNCellsInSupMod(geom.fNCellsInSupMod),
148     fNCellsInTower(geom.fNCellsInTower),
149     fNTRU(geom.fNTRU),
150     fNTRUEta(geom.fNTRUEta),
151     fNTRUPhi(geom.fNTRUPhi),
152     fTrd1Angle(geom.fTrd1Angle),
153     f2Trd1Dx2(geom.f2Trd1Dx2),
154     fPhiGapForSM(geom.fPhiGapForSM),
155     fKey110DEG(geom.fKey110DEG),
156     fPhiBoundariesOfSM(geom.fPhiBoundariesOfSM),
157     fPhiCentersOfSM(geom.fPhiCentersOfSM),
158     fEtaMaxOfTRD1(geom.fEtaMaxOfTRD1),
159     fTrd2AngleY(geom.fTrd2AngleY),
160     f2Trd2Dy2(geom.f2Trd2Dy2),
161     fEmptySpace(geom.fEmptySpace),
162     fTubsR(geom.fTubsR),
163     fTubsTurnAngle(geom.fTubsTurnAngle),
164     fCentersOfCellsEtaDir(geom.fCentersOfCellsEtaDir),
165     fCentersOfCellsXDir(geom.fCentersOfCellsXDir),
166     fCentersOfCellsPhiDir(geom.fCentersOfCellsPhiDir),
167     fEtaCentersOfCells(geom.fEtaCentersOfCells),
168     fPhiCentersOfCells(geom.fPhiCentersOfCells),
169     fShishKebabTrd1Modules(geom.fShishKebabTrd1Modules),
170     fNAdditionalOpts(geom.fNAdditionalOpts)
171 {
172   //copy ctor
173 }
174
175 //______________________________________________________________________
176 AliEMCALGeometry::~AliEMCALGeometry(void){
177     // dtor
178 }
179 //______________________________________________________________________
180 void AliEMCALGeometry::Init(void){
181   // Initializes the EMCAL parameters
182   // naming convention : GUV_WX_N_ gives the composition of a tower
183   // WX inform about the composition of the EM calorimeter section: 
184   //   thickness in mm of Pb radiator (W) and of scintillator (X), and number of scintillator layers (N)
185   // New geometry: EMCAL_55_25
186   // 24-aug-04 for shish-kebab
187   // SHISH_25 or SHISH_62
188   // 11-oct-05   - correction for pre final design
189   // Feb 06,2006 - decrease the weight of EMCAL
190   //
191   // Oct 30,2006 - SHISH_TRD1_CURRENT_1X1, SHISH_TRD1_CURRENT_2X2 or SHISH_TRD1_CURRENT_3X3;
192   //
193
194   fAdditionalOpts[0] = "nl=";    // number of sampling layers (fNECLayers)
195   fAdditionalOpts[1] = "pbTh=";  // cm, Thickness of the Pb   (fECPbRadThick)
196   fAdditionalOpts[2] = "scTh=";  // cm, Thickness of the Sc    (fECScintThick)
197   fAdditionalOpts[3] = "latSS=";  // cm, Thickness of lateral steel strip (fLateralSteelStrip)
198
199   fNAdditionalOpts = sizeof(fAdditionalOpts) / sizeof(char*);
200
201   fgInit = kFALSE; // Assume failed until proven otherwise.
202   fGeoName   = GetName();
203   fGeoName.ToUpper();
204   fKey110DEG = 0;
205   if(fGeoName.Contains("110DEG") || fGeoName.Contains("CURRENT")) fKey110DEG = 1; // for GetAbsCellId
206   fShishKebabTrd1Modules = 0;
207   fTrd2AngleY = f2Trd2Dy2 = fEmptySpace = fTubsR = fTubsTurnAngle = 0;
208
209   fNZ             = 114;        // granularity along Z (eta) 
210   fNPhi           = 168;        // granularity in phi (azimuth)
211   fArm1PhiMin     = 80.0;       // degrees, Starting EMCAL Phi position
212   fArm1PhiMax     = 190.0;      // degrees, Ending EMCAL Phi position
213   fArm1EtaMin     = -0.7;       // pseudorapidity, Starting EMCAL Eta position
214   fArm1EtaMax     = +0.7;       // pseudorapidity, Ending EMCAL Eta position
215   fIPDistance     = 454.0;      // cm, Radial distance to inner surface of EMCAL
216   fPhiGapForSM    = 0.;         // cm, only for final TRD1 geometry
217   for(int i=0; i<12; i++) fMatrixOfSM[i] = 0;
218
219   // geometry
220   if(fGeoName.Contains("SHISH")){ // Only shahslyk now
221     // 7-sep-05; integration issue
222     fArm1PhiMin     = 80.0;     // 60  -> 80
223     fArm1PhiMax     = 180.0;    // 180 -> 190
224
225     fNumberOfSuperModules = 10; // 12 = 6 * 2 (6 in phi, 2 in Z);
226     fSteelFrontThick = 2.54;    //  9-sep-04
227     fIPDistance      = 460.0;
228     fFrontSteelStrip = fPassiveScintThick = 0.0; // 13-may-05
229     fLateralSteelStrip = 0.025; // before MAY 2005 
230     fPhiModuleSize   = fEtaModuleSize   = 11.4;
231     fPhiTileSize = fEtaTileSize      = 5.52; // (11.4-5.52*2)/2. = 0.18 cm (wall thickness)
232     fNPhi            = 14;
233     fNZ              = 30;
234     fAlFrontThick    = fGap2Active = 0;
235     fNPHIdiv = fNETAdiv = 2;
236
237     fNECLayers       = 62;
238     fECScintThick    = fECPbRadThickness = 0.2;
239     fSampling        = 1.;  // 30-aug-04 - should be calculated
240     if(fGeoName.Contains("TWIST")) { // all about EMCAL module
241       fNZ             = 27;  // 16-sep-04
242     } else if(fGeoName.Contains("TRD")) {
243       fIPDistance      = 428.0;  //  11-may-05
244       fSteelFrontThick = 0.0;    // 3.17 -> 0.0; 28-mar-05 : no stell plate
245       fNPhi            = 12;
246       fSampling       = 12.327;
247       fPhiModuleSize = fEtaModuleSize = 12.26;
248       fNZ            = 26;     // 11-oct-04
249       fTrd1Angle     = 1.3;    // in degree
250 // 18-nov-04; 1./0.08112=12.327
251 // http://pdsfweb01.nersc.gov/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
252       if(fGeoName.Contains("TRD1")) {       // 30-jan-05
253         // for final design
254         fPhiGapForSM    = 2.;         // cm, only for final TRD1 geometry
255         if(fGeoName.Contains("MAY05") || fGeoName.Contains("WSUC") || fGeoName.Contains("FINAL") || fGeoName.Contains("CURRENT")){
256           fNumberOfSuperModules = 12; // 20-may-05
257           if(fGeoName.Contains("WSUC")) fNumberOfSuperModules = 1; // 27-may-05
258           fNECLayers     = 77;       // (13-may-05 from V.Petrov)
259           fPhiModuleSize = 12.5;     // 20-may-05 - rectangular shape
260           fEtaModuleSize = 11.9;
261           fECScintThick  = fECPbRadThickness = 0.16;// (13-may-05 from V.Petrov)
262           fFrontSteelStrip   = 0.025;// 0.025cm = 0.25mm  (13-may-05 from V.Petrov)
263           fLateralSteelStrip = 0.01; // 0.01cm  = 0.1mm   (13-may-05 from V.Petrov) - was 0.025
264           fPassiveScintThick = 0.8;  // 0.8cm   = 8mm     (13-may-05 from V.Petrov)
265           fNZ                = 24;
266           fTrd1Angle         = 1.5;  // 1.3 or 1.5
267
268           if(fGeoName.Contains("FINAL") || fGeoName.Contains("CURRENT")) { // 9-sep-05
269             fNumberOfSuperModules = 10;
270             if(GetKey110DEG()) {
271               fNumberOfSuperModules = 12;// last two modules have size 10 degree in phi (180<phi<190)
272               fArm1PhiMax = 200.0;       // for XEN1 and turn angle of super modules
273             }
274             if(fGeoName.Contains("FINAL")) {
275               fPhiModuleSize = 12.26 - fPhiGapForSM / Float_t(fNPhi); // first assumption
276             } else if(fGeoName.Contains("CURRENT")) {
277               fECScintThick      = 0.176; // 10% of weight reduction
278               fECPbRadThickness  = 0.144; //
279               fLateralSteelStrip = 0.015; // 0.015cm  = 0.15mm (Oct 30, from Fred)
280               fPhiModuleSize     = 12.00;
281               fPhiGapForSM       = (12.26 - fPhiModuleSize)*fNPhi; // have to check
282             }
283             fEtaModuleSize = fPhiModuleSize;
284             if(fGeoName.Contains("HUGE")) fNECLayers *= 3; // 28-oct-05 for analysing leakage    
285           }
286         }
287       } else if(fGeoName.Contains("TRD2")) {       // 30-jan-05
288         fSteelFrontThick = 0.0;         // 11-mar-05
289         fIPDistance+= fSteelFrontThick; // 1-feb-05 - compensate absence of steel plate
290         fTrd1Angle  = 1.64;             // 1.3->1.64
291         fTrd2AngleY = fTrd1Angle;       //  symmetric case now
292         fEmptySpace    = 0.2; // 2 mm
293         fTubsR         = fIPDistance; // 31-jan-05 - as for Fred case
294
295         fPhiModuleSize  = fTubsR*2.*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.);
296         fPhiModuleSize -= fEmptySpace/2.; // 11-mar-05  
297         fEtaModuleSize  = fPhiModuleSize; // 20-may-05 
298         fTubsTurnAngle  = 3.;
299       }
300       fNPHIdiv = fNETAdiv  = 2;   // 13-oct-04 - division again
301       if(fGeoName.Contains("3X3")) {   // 23-nov-04
302         fNPHIdiv = fNETAdiv  = 3;
303       } else if(fGeoName.Contains("4X4")) {
304         fNPHIdiv = fNETAdiv  = 4;
305       }
306     }
307     if(fGeoName.Contains("25")){
308       fNECLayers     = 25;
309       fECScintThick  = fECPbRadThickness = 0.5;
310     }
311     if(fGeoName.Contains("WSUC")){ // 18-may-05 - about common structure
312       fShellThickness = 30.; // should be change 
313       fNPhi = fNZ = 4; 
314     }
315
316     CheckAdditionalOptions();
317     DefineSamplingFraction();
318
319     fPhiTileSize = fPhiModuleSize/double(fNPHIdiv) - fLateralSteelStrip; // 13-may-05 
320     fEtaTileSize = fEtaModuleSize/double(fNETAdiv) - fLateralSteelStrip; // 13-may-05 
321
322     // constant for transition absid <--> indexes
323     fNCellsInTower  = fNPHIdiv*fNETAdiv;
324     fNCellsInSupMod = fNCellsInTower*fNPhi*fNZ;
325     fNCells         = fNCellsInSupMod*fNumberOfSuperModules;
326     if(GetKey110DEG()) fNCells -= fNCellsInSupMod;
327
328     fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness);
329     if(fGeoName.Contains("MAY05")) fLongModuleSize += (fFrontSteelStrip + fPassiveScintThick);
330
331     // 30-sep-04
332     if(fGeoName.Contains("TRD")) {
333       f2Trd1Dx2 = fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.);
334       if(fGeoName.Contains("TRD2")) {  // 27-jan-05
335         f2Trd2Dy2 = fPhiModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.);
336       }
337     }
338   } else Fatal("Init", "%s is an undefined geometry!", fGeoName.Data()) ; 
339
340   fNPhiSuperModule = fNumberOfSuperModules/2;
341   if(fNPhiSuperModule<1) fNPhiSuperModule = 1;
342
343   fShellThickness = fAlFrontThick + fGap2Active + fNECLayers*GetECScintThick()+(fNECLayers-1)*GetECPbRadThick();
344   if(fGeoName.Contains("SHISH")) {
345     fShellThickness = fSteelFrontThick + fLongModuleSize;
346     if(fGeoName.Contains("TWIST")) { // 13-sep-04
347       fShellThickness  = TMath::Sqrt(fLongModuleSize*fLongModuleSize + fPhiModuleSize*fEtaModuleSize);
348       fShellThickness += fSteelFrontThick;
349     } else if(fGeoName.Contains("TRD")) { // 1-oct-04
350       fShellThickness  = TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2);
351       fShellThickness += fSteelFrontThick;
352       // Local coordinates
353       fParSM[0] = GetShellThickness()/2.;        
354       fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
355       fParSM[2] = 350./2.;
356     }
357   }
358
359   fZLength        = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax); // Z coverage
360   fEnvelop[0]     = fIPDistance; // mother volume inner radius
361   fEnvelop[1]     = fIPDistance + fShellThickness; // mother volume outer r.
362   fEnvelop[2]     = 1.00001*fZLength; // add some padding for mother volume. 
363
364   fNumberOfSuperModules = 12;
365
366   // SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006 
367   fPhiBoundariesOfSM.Set(fNumberOfSuperModules);
368   fPhiCentersOfSM.Set(fNumberOfSuperModules/2);
369   fPhiBoundariesOfSM[0] = TMath::PiOver2() - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules)
370   fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance);
371   fPhiCentersOfSM[0]     = TMath::PiOver2();
372   for(int i=1; i<=4; i++) { // from 2th ro 9th
373     fPhiBoundariesOfSM[2*i]   = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i;
374     fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i;
375     fPhiCentersOfSM[i]         = fPhiCentersOfSM[0]     + 20.*TMath::DegToRad()*i;
376   }
377   fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad();
378   fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance);
379   fPhiCentersOfSM[5]      = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.; 
380
381   fgInit = kTRUE; 
382   
383   //TRU parameters. These parameters values are not the final ones.
384   fNTRU    = 3 ;
385   fNTRUEta = 3 ;
386   fNTRUPhi = 1 ;
387
388 }
389
390 void AliEMCALGeometry::PrintGeometry()
391 {
392   // Separate routine is callable from broswer; Nov 7,2006
393   printf("\nInit: geometry of EMCAL named %s is as follows:\n", fGeoName.Data());
394   printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
395   printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f -> for EMCAL envelope only\n",  
396            GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() );
397
398   printf( "               ECAL      : %d x (%f cm Pb, %f cm Sc) \n", 
399   GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ; 
400   printf("                fSampling %5.2f \n",  fSampling );
401   if(fGeoName.Contains("SHISH")){
402     printf(" fIPDistance       %6.3f cm \n", fIPDistance);
403     if(fSteelFrontThick>0.) 
404     printf(" fSteelFrontThick  %6.3f cm \n", fSteelFrontThick);
405     printf(" fNPhi %i   |  fNZ %i \n", fNPhi, fNZ);
406     printf(" fNCellsInTower %i : fNCellsInSupMod %i : fNCells %i\n",fNCellsInTower, fNCellsInSupMod, fNCells);
407     if(fGeoName.Contains("MAY05")){
408       printf(" fFrontSteelStrip         %6.4f cm (thickness of front steel strip)\n", 
409       fFrontSteelStrip);
410       printf(" fLateralSteelStrip       %6.4f cm (thickness of lateral steel strip)\n", 
411       fLateralSteelStrip);
412       printf(" fPassiveScintThick  %6.4f cm (thickness of front passive Sc tile)\n",
413       fPassiveScintThick);
414     }
415     printf(" X:Y module size     %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize);
416     printf(" X:Y   tile size     %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize);
417     printf(" #of sampling layers %i(fNECLayers) \n", fNECLayers);
418     printf(" fLongModuleSize     %6.3f cm \n", fLongModuleSize);
419     printf(" #supermodule in phi direction %i \n", fNPhiSuperModule );
420   }
421   if(fGeoName.Contains("TRD")) {
422     printf(" fTrd1Angle %7.4f\n", fTrd1Angle);
423     printf(" f2Trd1Dx2  %7.4f\n",  f2Trd1Dx2);
424     if(fGeoName.Contains("TRD2")) {
425       printf(" fTrd2AngleY     %7.4f\n", fTrd2AngleY);
426       printf(" f2Trd2Dy2       %7.4f\n", f2Trd2Dy2);
427       printf(" fTubsR          %7.2f cm\n", fTubsR);
428       printf(" fTubsTurnAngle  %7.4f\n", fTubsTurnAngle);
429       printf(" fEmptySpace     %7.4f cm\n", fEmptySpace);
430     } else if(fGeoName.Contains("TRD1")){
431       printf("SM dimensions(TRD1) : dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n", 
432       fParSM[0],fParSM[1],fParSM[2]);
433       printf(" fPhiGapForSM  %7.4f cm (%7.4f <- phi size in degree)\n",  
434       fPhiGapForSM, TMath::ATan2(fPhiGapForSM,fIPDistance)*TMath::RadToDeg());
435       if(GetKey110DEG()) printf(" Last two modules have size 10 degree in  phi (180<phi<190)\n");
436       printf(" phi SM boundaries \n"); 
437       for(int i=0; i<fPhiBoundariesOfSM.GetSize()/2.; i++) {
438         printf(" %i : %7.5f(%7.2f) -> %7.5f(%7.2f) : center %7.5f(%7.2f) \n", i, 
439         fPhiBoundariesOfSM[2*i], fPhiBoundariesOfSM[2*i]*TMath::RadToDeg(),
440                fPhiBoundariesOfSM[2*i+1], fPhiBoundariesOfSM[2*i+1]*TMath::RadToDeg(),
441                fPhiCentersOfSM[i], fPhiCentersOfSM[i]*TMath::RadToDeg());
442       }
443       printf(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n", 
444                fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1);
445
446       printf("\n Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize());
447       for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
448         printf(" ind %2.2i : z %8.3f : x %8.3f \n", i, 
449                fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i));
450         int ind=0; // Nov 21,2006
451         for(Int_t iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
452           ind = iphi*fCentersOfCellsEtaDir.GetSize() + i;
453           printf("%6.4f ", fEtaCentersOfCells[ind]);
454           if((iphi+1)%12 == 0) printf("\n");
455         }
456         printf("\n");
457       }
458
459       printf("\n Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize());
460       for(Int_t i=0; i<fCentersOfCellsPhiDir.GetSize(); i++) {
461         double phi=fPhiCentersOfCells.At(i);
462         printf(" ind %2.2i : y %8.3f : phi %7.5f(%6.2f) \n", i, fCentersOfCellsPhiDir.At(i), 
463                phi, phi*TMath::RadToDeg());
464       }
465     }
466   }
467 }
468
469 void AliEMCALGeometry::PrintCellIndexes(Int_t absId, int pri, char *tit)
470 {
471   // Service methods
472   Int_t nSupMod, nTower, nIphi, nIeta;
473   Int_t iphi, ieta;
474   TVector3 vg;
475
476   GetCellIndex(absId,  nSupMod, nTower, nIphi, nIeta);
477   printf(" %s | absId : %i -> nSupMod %i nTower %i nIphi %i nIeta %i \n", tit, absId,  nSupMod, nTower, nIphi, nIeta);
478   if(pri>0) {
479     GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
480     printf(" local SM index : iphi %i : ieta %i \n", iphi,ieta);
481     GetGlobal(absId, vg);
482     printf(" vglob : mag %7.2f : perp %7.2f : z %7.2f : eta %6.4f : phi %6.4f(%6.2f) \n", 
483            vg.Mag(), vg.Perp(), vg.Z(), vg.Eta(), vg.Phi(), vg.Phi()*TMath::RadToDeg());
484   }
485 }
486
487 //______________________________________________________________________
488 void AliEMCALGeometry::CheckAdditionalOptions()
489 {
490   // Feb 06,2006
491   //Additional options that
492   //can be used to select
493   //the specific geometry of 
494   //EMCAL to run
495
496   fArrayOpts = new TObjArray;
497   Int_t nopt = AliEMCALHistoUtilities::ParseString(fGeoName, *fArrayOpts);
498   if(nopt==1) { // no aditional option(s)
499     fArrayOpts->Delete();
500     delete fArrayOpts;
501     fArrayOpts = 0; 
502     return;
503   }              
504   for(Int_t i=1; i<nopt; i++){
505     TObjString *o = (TObjString*)fArrayOpts->At(i); 
506
507     TString addOpt = o->String();
508     Int_t indj=-1;
509     for(Int_t j=0; j<fNAdditionalOpts; j++) {
510       TString opt = fAdditionalOpts[j];
511       if(addOpt.Contains(opt,TString::kIgnoreCase)) {
512           indj = j;
513         break;
514       }
515     }
516     if(indj<0) {
517       AliDebug(2,Form("<E> option |%s| unavailable : ** look to the file AliEMCALGeometry.h **\n", 
518                       addOpt.Data()));
519       assert(0);
520     } else {
521       AliDebug(2,Form("<I> option |%s| is valid : number %i : |%s|\n", 
522                       addOpt.Data(), indj, fAdditionalOpts[indj]));
523       if       (addOpt.Contains("NL=",TString::kIgnoreCase))   {// number of sampling layers
524         sscanf(addOpt.Data(),"NL=%i", &fNECLayers);
525         AliDebug(2,Form(" fNECLayers %i (new) \n", fNECLayers));
526       } else if(addOpt.Contains("PBTH=",TString::kIgnoreCase)) {//Thickness of the Pb(fECPbRadThicknes)
527         sscanf(addOpt.Data(),"PBTH=%f", &fECPbRadThickness);
528       } else if(addOpt.Contains("SCTH=",TString::kIgnoreCase)) {//Thickness of the Sc(fECScintThick)
529         sscanf(addOpt.Data(),"SCTH=%f", &fECScintThick);
530       } else if(addOpt.Contains("LATSS=",TString::kIgnoreCase)) {// Thickness of lateral steel strip (fLateralSteelStrip)
531         sscanf(addOpt.Data(),"LATSS=%f", &fLateralSteelStrip);
532         AliDebug(2,Form(" fLateralSteelStrip %f (new) \n", fLateralSteelStrip));
533       }
534     }
535   }
536 }
537
538 void AliEMCALGeometry::DefineSamplingFraction()
539 {
540   // Jun 05,2006
541   // Look http://rhic.physics.wayne.edu/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
542   // Keep for compatibilty
543   //
544   if(fNECLayers == 69) {        // 10% layer reduction
545     fSampling = 12.55;
546   } else if(fNECLayers == 61) { // 20% layer reduction
547     fSampling = 12.80;
548   } else if(fNECLayers == 77) {
549     if       (fECScintThick>0.175 && fECScintThick<0.177) { // 10% Pb thicknes reduction
550       fSampling = 10.5; // fECScintThick = 0.176, fECPbRadThickness=0.144;
551     } else if(fECScintThick>0.191 && fECScintThick<0.193) { // 20% Pb thicknes reduction
552       fSampling = 8.93; // fECScintThick = 0.192, fECPbRadThickness=0.128;
553     }
554   }
555 }
556
557 //____________________________________________________________________________
558 void AliEMCALGeometry::FillTRU(const TClonesArray * digits, TClonesArray * ampmatrix, TClonesArray * timeRmatrix) {
559
560
561 //  Orders digits ampitudes list in fNTRU TRUs (384 cells) per supermodule. 
562 //  Each TRU is a TMatrixD, and they are kept in TClonesArrays. The number of 
563 //  TRU in phi is fNTRUPhi, and the number of TRU in eta is fNTRUEta.
564 //  Last 2 modules are half size in Phi, I considered that the number of TRU
565 //  is maintained for the last modules but decision not taken. If different, 
566 //  then this must be changed. 
567  
568
569   //Check data members
570
571   if(fNTRUEta*fNTRUPhi != fNTRU)
572     Error("FillTRU"," Wrong number of TRUS per Eta or Phi");
573
574   //Initilize and declare variables
575   //List of TRU matrices initialized to 0.
576   Int_t nCellsPhi  = fNPhi*2/fNTRUPhi;
577   Int_t nCellsPhi2 = fNPhi/fNTRUPhi; //HalfSize modules
578   Int_t nCellsEta  = fNZ*2/fNTRUEta;
579   Int_t id      = -1; 
580   Float_t amp   = -1;
581   Float_t timeR = -1;
582   Int_t iSupMod = -1;
583   Int_t nTower  = -1;
584   Int_t nIphi   = -1;
585   Int_t nIeta   = -1;
586   Int_t iphi    = -1;
587   Int_t ieta    = -1;
588
589   //List of TRU matrices initialized to 0.
590   for(Int_t k = 0; k < fNTRU*fNumberOfSuperModules; k++){
591     TMatrixD  * amptrus   = new TMatrixD(nCellsPhi,nCellsEta) ;
592     TMatrixD  * timeRtrus = new TMatrixD(nCellsPhi,nCellsEta) ;
593     for(Int_t i = 0; i < nCellsPhi; i++){
594       for(Int_t j = 0; j < nCellsEta; j++){
595         (*amptrus)(i,j) = 0.0;
596         (*timeRtrus)(i,j) = 0.0;
597       }
598     }
599     new((*ampmatrix)[k])   TMatrixD(*amptrus) ;
600     new((*timeRmatrix)[k]) TMatrixD(*timeRtrus) ; 
601   }
602   
603   AliEMCALDigit * dig ;
604   
605   //Digits loop to fill TRU matrices with amplitudes.
606   for(Int_t idig = 0 ; idig < digits->GetEntriesFast() ; idig++){
607     
608     dig = dynamic_cast<AliEMCALDigit *>(digits->At(idig)) ;
609     amp    = dig->GetAmp() ;   // Energy of the digit (arbitrary units)
610     id     = dig->GetId() ;    // Id label of the cell
611     timeR  = dig->GetTimeR() ; // Earliest time of the digit
612    
613     //Get eta and phi cell position in supermodule
614     Bool_t bCell = GetCellIndex(id, iSupMod, nTower, nIphi, nIeta) ;
615     if(!bCell)
616       Error("FillTRU","Wrong cell id number") ;
617     
618     GetCellPhiEtaIndexInSModule(iSupMod,nTower,nIphi, nIeta,iphi,ieta);
619
620     //Check to which TRU in the supermodule belongs the cell. 
621     //Supermodules are divided in a TRU matrix of dimension 
622     //(fNTRUPhi,fNTRUEta).
623     //Each TRU is a cell matrix of dimension (nCellsPhi,nCellsEta)
624
625     //First calculate the row and column in the supermodule 
626     //of the TRU to which the cell belongs.
627     Int_t col   = ieta/nCellsEta; 
628     Int_t row   = iphi/nCellsPhi; 
629     if(iSupMod > 9)
630       row   = iphi/nCellsPhi2; 
631     //Calculate label number of the TRU
632     Int_t itru  = row + col*fNTRUPhi + iSupMod*fNTRU ;  
633  
634     //Fill TRU matrix with cell values
635     TMatrixD * amptrus   = dynamic_cast<TMatrixD *>(ampmatrix->At(itru)) ;
636     TMatrixD * timeRtrus = dynamic_cast<TMatrixD *>(timeRmatrix->At(itru)) ;
637
638     //Calculate row and column of the cell inside the TRU with number itru
639     Int_t irow = iphi - row *  nCellsPhi;
640     if(iSupMod > 9)
641       irow = iphi - row *  nCellsPhi2;
642     Int_t icol = ieta - col *  nCellsEta;
643     
644     (*amptrus)(irow,icol) = amp ;
645     (*timeRtrus)(irow,icol) = timeR ;
646
647   }
648 }
649
650 //______________________________________________________________________
651 void AliEMCALGeometry::GetCellPhiEtaIndexInSModuleFromTRUIndex(const Int_t itru, const Int_t iphitru, const Int_t ietatru, Int_t &iphiSM, Int_t &ietaSM) const 
652 {
653   
654   // This method transforms the (eta,phi) index of cells in a 
655   // TRU matrix into Super Module (eta,phi) index.
656   
657   // Calculate in which row and column where the TRU are 
658   // ordered in the SM
659
660   Int_t col = itru/ fNTRUPhi ;
661   Int_t row = itru - col*fNTRUPhi ;
662    
663   //Calculate the (eta,phi) index in SM
664   Int_t nCellsPhi = fNPhi*2/fNTRUPhi;
665   Int_t nCellsEta = fNZ*2/fNTRUEta;
666   
667   iphiSM = nCellsPhi*row + iphitru  ;
668   ietaSM = nCellsEta*col + ietatru  ; 
669 }
670
671 //______________________________________________________________________
672 AliEMCALGeometry *  AliEMCALGeometry::GetInstance(){ 
673   // Returns the pointer of the unique instance
674   
675   AliEMCALGeometry * rv = static_cast<AliEMCALGeometry *>( fgGeom );
676   return rv; 
677 }
678
679 //______________________________________________________________________
680 AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
681                                                 const Text_t* title){
682     // Returns the pointer of the unique instance
683
684     AliEMCALGeometry * rv = 0; 
685     if ( fgGeom == 0 ) {
686         if ( strcmp(name,"") == 0 ) rv = 0;
687         else {
688             fgGeom = new AliEMCALGeometry(name, title);
689             if ( fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
690             else {
691                 rv = 0; 
692                 delete fgGeom; 
693                 fgGeom = 0; 
694             } // end if fgInit
695         } // end if strcmp(name,"")
696     }else{
697         if ( strcmp(fgGeom->GetName(), name) != 0) {
698           printf("\ncurrent geometry is %s : ", fgGeom->GetName());
699           printf(" you cannot call %s ", name);  
700         }else{
701           rv = (AliEMCALGeometry *) fgGeom; 
702         } // end 
703     }  // end if fgGeom
704     return rv; 
705 }
706
707 Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const {
708   // Checks whether point is inside the EMCal volume, used in AliEMCALv*.cxx
709   //
710   // Code uses cylindrical approximation made of inner radius (for speed)
711   //
712   // Points behind EMCAl, i.e. R > outer radius, but eta, phi in acceptance 
713   // are considered to inside
714
715   Double_t r=sqrt(x*x+y*y);
716
717   if ( r > fEnvelop[0] ) {
718      Double_t theta;
719      theta  =    TMath::ATan2(r,z);
720      Double_t eta;
721      if(theta == 0) 
722        eta = 9999;
723      else 
724        eta    =   -TMath::Log(TMath::Tan(theta/2.));
725      if (eta < fArm1EtaMin || eta > fArm1EtaMax)
726        return 0;
727  
728      Double_t phi = TMath::ATan2(y,x) * 180./TMath::Pi();
729      if (phi > fArm1PhiMin && phi < fArm1PhiMax)
730        return 1;
731   }
732   return 0;
733 }
734 // ==
735
736 //
737 // == Shish-kebab cases ==
738 //
739 Int_t AliEMCALGeometry::GetAbsCellId(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta) const
740
741   // 27-aug-04; 
742   // corr. 21-sep-04; 
743   //       13-oct-05; 110 degree case
744   // May 31, 2006; ALICE numbering scheme:
745   // 0 <= nSupMod < fNumberOfSuperModules
746   // 0 <= nTower  < fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1)
747   // 0 <= nIphi   < fNPHIdiv
748   // 0 <= nIeta   < fNETAdiv
749   // 0 <= absid   < fNCells
750   static Int_t id=0; // have to change from 0 to fNCells-1
751   if(fKey110DEG == 1 && nSupMod >= 10) { // 110 degree case; last two supermodules
752     id  = fNCellsInSupMod*10 + (fNCellsInSupMod/2)*(nSupMod-10);
753   } else {
754     id  = fNCellsInSupMod*nSupMod;
755   }
756   id += fNCellsInTower *nTower;
757   id += fNPHIdiv *nIphi;
758   id += nIeta;
759   if(id<0 || id >= fNCells) {
760 //     printf(" wrong numerations !!\n");
761 //     printf("    id      %6i(will be force to -1)\n", id);
762 //     printf("    fNCells %6i\n", fNCells);
763 //     printf("    nSupMod %6i\n", nSupMod);
764 //     printf("    nTower  %6i\n", nTower);
765 //     printf("    nIphi   %6i\n", nIphi);
766 //     printf("    nIeta   %6i\n", nIeta);
767     id = -TMath::Abs(id); // if negative something wrong
768   }
769   return id;
770 }
771
772 Bool_t  AliEMCALGeometry::CheckAbsCellId(Int_t absId) const
773
774   // May 31, 2006; only trd1 now
775   if(absId<0 || absId >= fNCells) return kFALSE;
776   else                            return kTRUE;
777 }
778
779 Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nTower,Int_t &nIphi,Int_t &nIeta) const
780
781   // 21-sep-04; 19-oct-05;
782   // May 31, 2006; ALICE numbering scheme:
783   static Int_t tmp=0, sm10=0;
784   if(!CheckAbsCellId(absId)) return kFALSE;
785
786   sm10 = fNCellsInSupMod*10;
787   if(fKey110DEG == 1 && absId >= sm10) { // 110 degree case; last two supermodules  
788     nSupMod = (absId-sm10) / (fNCellsInSupMod/2) + 10;
789     tmp     = (absId-sm10) % (fNCellsInSupMod/2);
790   } else {
791     nSupMod = absId / fNCellsInSupMod;
792     tmp     = absId % fNCellsInSupMod;
793   }
794
795   nTower  = tmp / fNCellsInTower;
796   tmp     = tmp % fNCellsInTower;
797   nIphi   = tmp / fNPHIdiv;
798   nIeta   = tmp % fNPHIdiv;
799
800   return kTRUE;
801 }
802
803 void AliEMCALGeometry::GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower,  int &iphim, int &ietam) const
804
805   // added nSupMod; - 19-oct-05 !
806   // Alice numbering scheme        - Jun 01,2006 
807   // ietam, iphi - indexes of module in two dimensional grid of SM
808   // ietam - have to change from 0 to fNZ-1
809   // iphim - have to change from 0 to nphi-1 (fNPhi-1 or fNPhi/2-1)
810   static Int_t nphi;
811
812   if(fKey110DEG == 1 && nSupMod>=10) nphi = fNPhi/2;
813   else                               nphi = fNPhi;
814
815   ietam = nTower/nphi;
816   iphim = nTower%nphi;
817 }
818
819 void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta, 
820 int &iphi, int &ieta) const
821
822   // 
823   // Added nSupMod; Nov 25, 05
824   // Alice numbering scheme  - Jun 01,2006 
825   // ieta, iphi - indexes of cell(tower) in two dimensional grid of SM
826   // ieta - have to change from 0 to (fNZ*fNETAdiv-1)
827   // iphi - have to change from 0 to (fNPhi*fNPHIdiv-1 or fNPhi*fNPHIdiv/2-1)
828   //
829   static Int_t iphim, ietam;
830
831   GetModulePhiEtaIndexInSModule(nSupMod,nTower, iphim, ietam); 
832   //  ieta  = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM) 
833   ieta  = ietam*fNETAdiv + (fNETAdiv - 1 - nIeta); // x(module) = -z(SM) 
834   iphi  = iphim*fNPHIdiv + nIphi;     // y(module) =  y(SM) 
835
836   if(iphi<0 || ieta<0)
837   AliDebug(1,Form(" nSupMod %i nTower %i nIphi %i nIeta %i => ieta %i iphi %i\n", 
838   nSupMod, nTower, nIphi, nIeta, ieta, iphi));
839 }
840
841 Int_t  AliEMCALGeometry::GetSuperModuleNumber(Int_t absId)  const
842 {
843   // Return the number of the  supermodule given the absolute
844   // ALICE numbering id
845
846   static Int_t nSupMod, nTower, nIphi, nIeta;
847   GetCellIndex(absId, nSupMod, nTower, nIphi, nIeta);
848   return nSupMod;
849
850
851 void  AliEMCALGeometry::GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta, 
852                         Int_t &iphim, Int_t &ietam, Int_t &nTower) const
853 {
854   // Transition from cell indexes (ieta,iphi) to module indexes (ietam,iphim, nTower)
855   static Int_t nphi;
856   nphi  = GetNumberOfModuleInPhiDirection(nSupMod);  
857
858   ietam  = ieta/fNETAdiv;
859   iphim  = iphi/fNPHIdiv;
860   nTower = ietam * nphi + iphim; 
861 }
862
863 Int_t  AliEMCALGeometry::GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const
864 {
865   // Transition from super module number(nSupMod) and cell indexes (ieta,iphi) to absId
866   static Int_t ietam, iphim, nTower;
867   static Int_t nIeta, nIphi; // cell indexes in module
868
869   GetModuleIndexesFromCellIndexesInSModule(nSupMod, iphi, ieta, ietam, iphim, nTower);
870
871   nIeta = ieta%fNETAdiv;
872   nIeta = fNETAdiv - 1 - nIeta;
873   nIphi = iphi%fNPHIdiv;
874
875   return GetAbsCellId(nSupMod, nTower, nIphi, nIeta);
876 }
877
878
879 // Methods for AliEMCALRecPoint - Feb 19, 2006
880 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const
881 {
882   // Look to see what the relative
883   // position inside a given cell is
884   // for a recpoint.
885   // Alice numbering scheme - Jun 08, 2006
886
887   static Int_t nSupMod, nTower, nIphi, nIeta, iphi, ieta;
888   static Int_t phiIndexShift=6;
889   if(!CheckAbsCellId(absId)) return kFALSE;
890
891   GetCellIndex(absId, nSupMod, nTower, nIphi, nIeta);
892   GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi, ieta); 
893  
894   xr = fCentersOfCellsXDir.At(ieta);
895   zr = fCentersOfCellsEtaDir.At(ieta);
896
897   if(nSupMod<10) {
898     yr = fCentersOfCellsPhiDir.At(iphi);
899   } else {
900     yr = fCentersOfCellsPhiDir.At(iphi + phiIndexShift);
901   }
902   if(absId==1104 || absId==1105) 
903   AliDebug(2,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
904
905   return kTRUE;
906 }
907
908 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t loc[3]) const
909 {
910   // Alice numbering scheme - Jun 03, 2006
911   loc[0] = loc[1] = loc[2]=0.0;
912   if(RelPosCellInSModule(absId, loc[0],loc[1],loc[2])) {
913     return kTRUE;
914   }
915   return kFALSE;
916 }
917
918 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, TVector3 &vloc) const
919 {
920   static Double_t loc[3];
921   if(RelPosCellInSModule(absId,loc)) {
922     vloc.SetXYZ(loc[0], loc[1], loc[2]);
923     return kTRUE;
924   } else {
925     vloc.SetXYZ(0,0,0);
926     return kFALSE;
927   }
928   // Alice numbering scheme - Jun 03, 2006
929 }
930
931 void AliEMCALGeometry::CreateListOfTrd1Modules()
932 {
933   // Generate the list of Trd1 modules
934   // which will make up the EMCAL
935   // geometry
936
937   AliDebug(2,Form(" AliEMCALGeometry::CreateListOfTrd1Modules() started "));
938
939   AliEMCALShishKebabTrd1Module *mod=0, *mTmp=0; // current module
940   if(fShishKebabTrd1Modules == 0) {
941     fShishKebabTrd1Modules = new TList;
942     fShishKebabTrd1Modules->SetName("ListOfTRD1");
943     for(int iz=0; iz< GetNZ(); iz++) { 
944       if(iz==0) { 
945         mod  = new AliEMCALShishKebabTrd1Module(TMath::Pi()/2.,this);
946       } else {
947         mTmp  = new AliEMCALShishKebabTrd1Module(*mod);
948         mod   = mTmp;
949       }
950       fShishKebabTrd1Modules->Add(mod);
951     }
952   } else {
953     AliDebug(2,Form(" Already exits : "));
954   }
955   mod = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(fShishKebabTrd1Modules->GetSize()-1);
956   fEtaMaxOfTRD1 = mod->GetMaxEtaOfModule(0);
957
958   AliDebug(2,Form(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n", 
959                   fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1));
960   // Feb 20,2006;
961   // Jun 01, 2006 - ALICE numbering scheme
962   // define grid for cells in eta(z) and x directions in local coordinates system of SM
963   // Works just for 2x2 case only -- ?? start here
964   // 
965   //
966   // Define grid for cells in phi(y) direction in local coordinates system of SM
967   // as for 2X2 as for 3X3 - Nov 8,2006
968   // 
969   AliDebug(2,Form(" Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize()));
970   Int_t ind=0; // this is phi index
971   Int_t iphi=0, ieta=0, nTower=0;
972   Double_t xr, zr, theta, phi, eta, r, x,y;
973   TVector3 vglob;
974   Double_t ytCenterModule, ytCenterCell;
975
976   fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
977   fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
978
979   Double_t R0 = GetIPDistance() + GetLongModuleSize()/2.;
980   for(Int_t it=0; it<fNPhi; it++) { // cycle on modules
981     ytCenterModule = -fParSM[1] + fPhiModuleSize*(2*it+1)/2;  // center of module
982     for(Int_t ic=0; ic<fNPHIdiv; ic++) { // cycle on cells in module
983       if(fNPHIdiv==2) {
984         ytCenterCell = ytCenterModule + fPhiTileSize *(2*ic-1)/2.;
985       } else if(fNPHIdiv==3){
986         ytCenterCell = ytCenterModule + fPhiTileSize *(ic-1);
987       }
988       fCentersOfCellsPhiDir.AddAt(ytCenterCell,ind);
989       // Define grid on phi direction
990       // Grid is not the same for different eta bin;
991       // Effect is small but is still here
992       phi = TMath::ATan2(ytCenterCell, R0);
993       fPhiCentersOfCells.AddAt(phi, ind);
994
995       AliDebug(2,Form(" ind %2.2i : y %8.3f ", ind, fCentersOfCellsPhiDir.At(ind))); 
996       ind++;
997     }
998   }
999
1000   fCentersOfCellsEtaDir.Set(fNZ *fNETAdiv);
1001   fCentersOfCellsXDir.Set(fNZ *fNETAdiv);
1002   fEtaCentersOfCells.Set(fNZ *fNETAdiv * fNPhi*fNPHIdiv);
1003   AliDebug(2,Form(" Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize()));
1004   for(Int_t it=0; it<fNZ; it++) {
1005     AliEMCALShishKebabTrd1Module *trd1 = GetShishKebabModule(it);
1006     nTower = fNPhi*it;
1007     for(Int_t ic=0; ic<fNETAdiv; ic++) {
1008       if(fNPHIdiv==2) {
1009         trd1->GetCenterOfCellInLocalCoordinateofSM(ic, xr, zr);    // case of 2X2
1010         GetCellPhiEtaIndexInSModule(0, nTower, 0, ic, iphi, ieta); // don't depend from phi - ieta in action
1011         fCentersOfCellsXDir.AddAt(float(xr) - fParSM[0],ieta);
1012         fCentersOfCellsEtaDir.AddAt(float(zr) - fParSM[2],ieta);
1013       } if(fNPHIdiv==3) {
1014         trd1->GetCenterOfCellInLocalCoordinateofSM_3X3(ic, xr, zr);  // case of 3X3
1015         GetCellPhiEtaIndexInSModule(0, nTower, 0, ic, iphi, ieta);   // don't depend from phi - ieta in action
1016         fCentersOfCellsXDir.AddAt(float(xr) - fParSM[0],ieta);
1017         fCentersOfCellsEtaDir.AddAt(float(zr) - fParSM[2],ieta);
1018       }
1019       // Define grid on eta direction for each bin in phi
1020       for(int iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
1021         x = xr + trd1->GetRadius();
1022         y = fCentersOfCellsPhiDir[iphi];
1023         r = TMath::Sqrt(x*x + y*y + zr*zr);
1024         theta = TMath::ACos(zr/r);
1025         eta   = AliEMCALShishKebabTrd1Module::ThetaToEta(theta);
1026         //        ind   = ieta*fCentersOfCellsPhiDir.GetSize() + iphi;
1027         ind   = iphi*fCentersOfCellsEtaDir.GetSize() + ieta;
1028         fEtaCentersOfCells.AddAt(eta, ind);
1029       }
1030       //printf(" ieta %i : xr + trd1->GetRadius() %f : zr %f : eta %f \n", ieta, xr + trd1->GetRadius(), zr, eta);
1031     }
1032   }
1033   for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1034     AliDebug(2,Form(" ind %2.2i : z %8.3f : x %8.3f", i+1, 
1035                     fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i)));
1036   }
1037
1038 }
1039
1040 void  AliEMCALGeometry::GetTransformationForSM()
1041 {
1042   //Uses the geometry manager to
1043   //load the transformation matrix
1044   //for the supermodules
1045
1046   static Bool_t transInit=kFALSE;
1047   if(transInit) return;
1048
1049   int i=0;
1050   if(gGeoManager == 0) {
1051     Info("CreateTransformationForSM() "," Load geometry : TGeoManager::Import()");
1052     assert(0);
1053   }
1054   TGeoNode *tn = gGeoManager->GetTopNode();
1055   TGeoNode *node=0, *xen1 = 0;
1056   for(i=0; i<tn->GetNdaughters(); i++) {
1057     node = tn->GetDaughter(i);
1058     TString ns(node->GetName());
1059     if(ns.Contains(GetNameOfEMCALEnvelope())) {
1060       xen1 = node;
1061       break;
1062     }
1063   }
1064   if(!xen1) {
1065     Info("CreateTransformationForSM() "," geometry has not EMCAL envelope with name %s", 
1066     GetNameOfEMCALEnvelope());
1067     assert(0);
1068   }
1069   printf(" i %i : EMCAL Envelope is %s : #SM %i \n", i, xen1->GetName(), xen1->GetNdaughters());
1070   for(i=0; i<xen1->GetNdaughters(); i++) {
1071     TGeoNodeMatrix *sm = (TGeoNodeMatrix*)xen1->GetDaughter(i);
1072     fMatrixOfSM[i] = sm->GetMatrix();
1073     //Compiler doesn't like this syntax...
1074     //    printf(" %i : matrix %x \n", i, fMatrixOfSM[i]);
1075   }
1076   transInit = kTRUE;
1077 }
1078
1079 void AliEMCALGeometry::GetGlobal(const Double_t *loc, Double_t *glob, int ind) const
1080 {
1081   // Figure out the global numbering
1082   // of a given supermodule from the
1083   // local numbering
1084   // Alice numbering - Jun 03,2006
1085   //  if(fMatrixOfSM[0] == 0) GetTransformationForSM();
1086
1087   if(ind>=0 && ind < GetNumberOfSuperModules()) {
1088     fMatrixOfSM[ind]->LocalToMaster(loc, glob);
1089   }
1090 }
1091
1092 void AliEMCALGeometry::GetGlobal(const TVector3 &vloc, TVector3 &vglob, int ind) const
1093 {
1094   //Figure out the global numbering
1095   //of a given supermodule from the
1096   //local numbering given a 3-vector location
1097
1098   static Double_t tglob[3], tloc[3];
1099   vloc.GetXYZ(tloc);
1100   GetGlobal(tloc, tglob, ind);
1101   vglob.SetXYZ(tglob[0], tglob[1], tglob[2]);
1102 }
1103
1104 void AliEMCALGeometry::GetGlobal(Int_t absId , double glob[3]) const
1105
1106   // Alice numbering scheme - Jun 03, 2006
1107   static Int_t nSupMod, nModule, nIphi, nIeta;
1108   static double loc[3];
1109
1110   glob[0]=glob[1]=glob[2]=0.0; // bad case
1111   if(RelPosCellInSModule(absId, loc)) {
1112     GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1113     fMatrixOfSM[nSupMod]->LocalToMaster(loc, glob);
1114   }
1115 }
1116
1117 void AliEMCALGeometry::GetGlobal(Int_t absId , TVector3 &vglob) const
1118
1119   // Alice numbering scheme - Jun 03, 2006
1120   static Double_t glob[3];
1121
1122   GetGlobal(absId, glob);
1123   vglob.SetXYZ(glob[0], glob[1], glob[2]);
1124
1125 }
1126
1127 void AliEMCALGeometry::GetGlobal(const AliRecPoint *rp, TVector3 &vglob) const
1128 {
1129   // Figure out the global numbering
1130   // of a given supermodule from the
1131   // local numbering for RecPoints
1132
1133   static TVector3 vloc;
1134   static Int_t nSupMod, nModule, nIphi, nIeta;
1135
1136   AliRecPoint *rpTmp = (AliRecPoint*)rp; // const_cast ??
1137   if(!rpTmp) return;
1138   AliEMCALRecPoint *rpEmc = (AliEMCALRecPoint*)rpTmp;
1139
1140   GetCellIndex(rpEmc->GetAbsId(0), nSupMod, nModule, nIphi, nIeta);
1141   rpTmp->GetLocalPosition(vloc);
1142   GetGlobal(vloc, vglob, nSupMod);
1143 }
1144
1145 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Double_t &eta,Double_t &phi) const
1146 {
1147   // Nov 16, 2006- float to double
1148   // version for TRD1 only
1149   static TVector3 vglob;
1150   GetGlobal(absId, vglob);
1151   eta = vglob.Eta();
1152   phi = vglob.Phi();
1153 }
1154
1155 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Float_t &eta,Float_t &phi) const
1156 {
1157   // Nov 16,2006 - should be discard in future
1158   static TVector3 vglob;
1159   GetGlobal(absId, vglob);
1160   eta = float(vglob.Eta());
1161   phi = float(vglob.Phi());
1162 }
1163
1164 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSM(Int_t nSupMod, Double_t &phiMin, Double_t &phiMax) const
1165 {
1166   // 0<= nSupMod <=11; phi in rad
1167   static int i;
1168   if(nSupMod<0 || nSupMod >11) return kFALSE; 
1169   i = nSupMod/2;
1170   phiMin = fPhiBoundariesOfSM[2*i];
1171   phiMax = fPhiBoundariesOfSM[2*i+1];
1172   return kTRUE; 
1173 }
1174
1175 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const
1176 {
1177   // 0<= nPhiSec <=4; phi in rad
1178   // 0;  gap boundaries between  0th&2th  | 1th&3th SM
1179   // 1;  gap boundaries between  2th&4th  | 3th&5th SM
1180   // 2;  gap boundaries between  4th&6th  | 5th&7th SM
1181   // 3;  gap boundaries between  6th&8th  | 7th&9th SM
1182   // 4;  gap boundaries between  8th&10th | 9th&11th SM
1183   if(nPhiSec<0 || nPhiSec >4) return kFALSE; 
1184   phiMin = fPhiBoundariesOfSM[2*nPhiSec+1];
1185   phiMax = fPhiBoundariesOfSM[2*nPhiSec+2];
1186   return kTRUE; 
1187 }
1188
1189 Bool_t AliEMCALGeometry::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const
1190
1191   // Return false if phi belongs a phi cracks between SM
1192  
1193   static Int_t i;
1194
1195   if(TMath::Abs(eta) > fEtaMaxOfTRD1) return kFALSE;
1196
1197   phi = TVector2::Phi_0_2pi(phi); // move phi to (0,2pi) boundaries
1198   for(i=0; i<6; i++) {
1199     if(phi>=fPhiBoundariesOfSM[2*i] && phi<=fPhiBoundariesOfSM[2*i+1]) {
1200       nSupMod = 2*i;
1201       if(eta < 0.0) nSupMod++;
1202       //      Info("SuperModuleNumberFromEtaPhi", Form(" nSupMod %i : eta %5.3f : phi %5.3f(%5.1f) ", 
1203       //                       nSupMod, eta, phi, phi*TMath::RadToDeg()));
1204       return kTRUE;
1205     }
1206   }
1207   AliDebug(1,Form("eta %f phi %f(%5.2f) : nSupMod %i : #bound %i", eta,phi,phi*TMath::RadToDeg(), nSupMod,i));
1208   return kFALSE;
1209 }
1210
1211 Bool_t AliEMCALGeometry::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi, Int_t &absId) const
1212 {
1213   // Nov 17,2006
1214   // stay here - phi problem as usual 
1215   static Int_t nSupMod, i, ieta, iphi, etaShift, nphi;
1216   static Double_t absEta=0.0, d=0.0, dmin=0.0, phiLoc;
1217   absId = nSupMod = - 1;
1218   if(SuperModuleNumberFromEtaPhi(eta, phi, nSupMod)) {
1219     // phi index first
1220     phi    = TVector2::Phi_0_2pi(phi);
1221     phiLoc = phi - fPhiCentersOfSM[nSupMod/2];
1222     nphi   = fPhiCentersOfCells.GetSize();
1223     if(nSupMod>=10) {
1224       phiLoc = phi - 190.*TMath::DegToRad();
1225       nphi  /= 2;
1226     }
1227
1228     dmin   = TMath::Abs(fPhiCentersOfCells[0]-phiLoc);
1229     iphi   = 0;
1230     for(i=1; i<nphi; i++) {
1231       d = TMath::Abs(fPhiCentersOfCells[i] - phiLoc);
1232       if(d < dmin) {
1233         dmin = d;
1234         iphi = i;
1235       }
1236       //      printf(" i %i : d %f : dmin %f : fPhiCentersOfCells[i] %f \n", i, d, dmin, fPhiCentersOfCells[i]);
1237     }
1238     // odd SM are turned with respect of even SM - reverse indexes
1239     AliDebug(2,Form(" iphi %i : dmin %f (phi %f, phiLoc %f ) ", iphi, dmin, phi, phiLoc));
1240     // eta index
1241     absEta   = TMath::Abs(eta);
1242     etaShift = iphi*fCentersOfCellsEtaDir.GetSize();
1243     dmin     = TMath::Abs(fEtaCentersOfCells[etaShift]-absEta);
1244     ieta     = 0;
1245     for(i=1; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1246       d = TMath::Abs(fEtaCentersOfCells[i+etaShift] - absEta);
1247       if(d < dmin) {
1248         dmin = d;
1249         ieta = i;
1250       }
1251     }
1252     AliDebug(2,Form(" ieta %i : dmin %f (eta=%f) : nSupMod %i ", ieta, dmin, eta, nSupMod));
1253
1254     if(eta<0) iphi = (nphi-1) - iphi;
1255     absId = GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);
1256
1257     return kTRUE;
1258   }
1259   return kFALSE;
1260 }
1261
1262 AliEMCALShishKebabTrd1Module* AliEMCALGeometry::GetShishKebabModule(Int_t neta)
1263 {
1264   //This method was too long to be
1265   //included in the header file - the
1266   //rule checker complained about it's
1267   //length, so we move it here.  It returns the
1268   //shishkebabmodule at a given eta index point.
1269
1270   static AliEMCALShishKebabTrd1Module* trd1=0;
1271   if(fShishKebabTrd1Modules && neta>=0 && neta<fShishKebabTrd1Modules->GetSize()) {
1272     trd1 = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(neta);
1273   } else trd1 = 0;
1274   return trd1;
1275 }
1276
1277 void AliEMCALGeometry::Browse(TBrowser* b)
1278 {
1279   if(fShishKebabTrd1Modules) b->Add(fShishKebabTrd1Modules);
1280 }
1281
1282 Bool_t AliEMCALGeometry::IsFolder() const
1283 {
1284   if(fShishKebabTrd1Modules) return kTRUE;
1285   else                       return kFALSE;
1286 }