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