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