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