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