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