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