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