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