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