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