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