1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
26 // EMCAL geometry tree:
27 // EMCAL -> superModule -> module -> tower(cell)
29 // absId -> nSupMod -> nModule -> (nIphi,nIeta)
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
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
45 //*-- Author: Sahal Yacoob (LBL / UCT)
46 // and : Yves Schutz (SUBATECH)
47 // and : Jennifer Klay (LBL)
48 // and : Aleksei Pavlinov (WSU)
53 // --- Root header files ---
54 #include <Riostream.h>
56 #include <TClonesArray.h>
57 #include <TGeoManager.h>
58 #include <TGeoMatrix.h>
62 #include <TObjArray.h>
63 #include <TObjString.h>
66 #include <TParticle.h>
71 #include "AliEMCALGeometry.h"
72 #include "AliEMCALShishKebabTrd1Module.h"
73 #include "AliEMCALRecPoint.h"
74 #include "AliEMCALDigit.h"
75 #include "AliEMCALHistoUtilities.h"
77 ClassImp(AliEMCALGeometry)
79 // these initialisations are needed for a singleton
80 AliEMCALGeometry *AliEMCALGeometry::fgGeom = 0;
81 Bool_t AliEMCALGeometry::fgInit = kFALSE;
82 Char_t* AliEMCALGeometry::fgDefaultGeometryName = "EMCAL_COMPLETE";
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.
89 // AliEMCALGeometry* g = AliEMCALGeometry::GetInstance(name,title); // first time
91 // g = AliEMCALGeometry::GetInstance(); // after first time
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");
99 AliEMCALGeometry::AliEMCALGeometry()
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),
110 fNModulesInTRUEta(0),
111 fNModulesInTRUPhi(0),
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.)
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
127 AliDebug(2, "AliEMCALGeometry : default ctor ");
129 //______________________________________________________________________
130 AliEMCALGeometry::AliEMCALGeometry(const Text_t* name, const Text_t* title)
131 : AliGeometry(name, title),
132 fGeoName(0),fArrayOpts(0),fNAdditionalOpts(0),fECPbRadThickness(0.),fECScintThick(0.),
133 fNECLayers(0),fArm1PhiMin(0.),fArm1PhiMax(0.),fArm1EtaMin(0.),fArm1EtaMax(0.),fIPDistance(0.),
134 fShellThickness(0.),fZLength(0.),fNZ(0),fNPhi(0),fSampling(0.),fNumberOfSuperModules(0),
135 fFrontSteelStrip(0.),fLateralSteelStrip(0.),fPassiveScintThick(0.),fPhiModuleSize(0.),
136 fEtaModuleSize(0.),fPhiTileSize(0.),fEtaTileSize(0.),fLongModuleSize(0.),fNPhiSuperModule(0),
137 fNPHIdiv(0),fNETAdiv(0), fNCells(0),fNCellsInSupMod(0),fNCellsInModule(0),
141 fNModulesInTRUEta(0),
142 fNModulesInTRUPhi(0),
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.)
154 // ctor only for internal usage (singleton)
155 AliDebug(2, Form("AliEMCALGeometry(%s,%s) ", name,title));
159 CreateListOfTrd1Modules();
161 if (AliDebugLevel()>=2) {
166 //______________________________________________________________________
167 AliEMCALGeometry::AliEMCALGeometry(const AliEMCALGeometry& 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),
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),
201 fNTRUEta(geom.fNTRUEta),
202 fNTRUPhi(geom.fNTRUPhi),
203 fNModulesInTRUEta(geom.fNModulesInTRUEta),
204 fNModulesInTRUPhi(geom.fNModulesInTRUPhi),
205 fNEtaSubOfTRU(geom.fNEtaSubOfTRU),
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),
229 fTubsTurnAngle(geom.fTubsTurnAngle)
234 //______________________________________________________________________
235 AliEMCALGeometry::~AliEMCALGeometry(void){
239 //______________________________________________________________________
240 void AliEMCALGeometry::Init(void){
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
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)
255 fNAdditionalOpts = sizeof(fAdditionalOpts) / sizeof(char*);
258 fgInit = kFALSE; // Assume failed until proven otherwise.
259 fGeoName = GetName();
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";
267 fGeoName = "EMCAL_PDC06";
270 if(fGeoName.Contains("WSUC")) fGeoName = "EMCAL_WSUC";
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()) ;
277 // Option to know whether we have the "half" supermodule(s) or not
279 if(fGeoName.Contains("COMPLETE") || fGeoName.Contains("PDC06")) fKey110DEG = 1; // for GetAbsCellId
280 fShishKebabTrd1Modules = 0;
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
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
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
307 fPhiModuleSize = 12.26 - fPhiGapForSM / Float_t(fNPhi); // first assumption
308 fEtaModuleSize = fPhiModuleSize;
310 fZLength = 700.; // Z coverage (cm)
313 //needs to be called for each geometry and before setting geometry
314 //parameters which can depend on the outcome
315 CheckAdditionalOptions();
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();
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
331 CheckAdditionalOptions();
334 // constant for transition absid <--> indexes
335 fNCellsInModule = fNPHIdiv*fNETAdiv;
336 fNCellsInSupMod = fNCellsInModule*fNPhi*fNZ;
337 fNCells = fNCellsInSupMod*fNumberOfSuperModules;
338 if(GetKey110DEG()) fNCells -= fNCellsInSupMod;
340 fNPhiSuperModule = fNumberOfSuperModules/2;
341 if(fNPhiSuperModule<1) fNPhiSuperModule = 1;
343 fPhiTileSize = fPhiModuleSize/double(fNPHIdiv) - fLateralSteelStrip; // 13-may-05
344 fEtaTileSize = fEtaModuleSize/double(fNETAdiv) - fLateralSteelStrip; // 13-may-05
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);
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
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
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;
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.;
381 //called after setting of scintillator and lead layer parameters
382 DefineSamplingFraction();
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
393 // 3*6*10 + 2*6*2 = 204 -> matrix (nphi(17), neta(12))
399 //___________________________________________________________________
400 void AliEMCALGeometry::PrintGeometry()
402 // Separate routine is callable from broswer; Nov 7,2006
403 printf("\nInit: geometry of EMCAL named %s :\n", fGeoName.Data());
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());
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() );
414 printf( " ECAL : %d x (%f cm Pb, %f cm Sc) \n",
415 GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ;
416 printf(" fSampling %5.2f \n", fSampling );
417 printf(" fIPDistance %6.3f cm \n", fIPDistance);
418 printf(" fNPhi %i | fNZ %i \n", fNPhi, fNZ);
419 printf(" fNCellsInModule %i : fNCellsInSupMod %i : fNCells %i\n",fNCellsInModule, fNCellsInSupMod, fNCells);
420 printf(" X:Y module size %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize);
421 printf(" X:Y tile size %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize);
422 printf(" #of sampling layers %i(fNECLayers) \n", fNECLayers);
423 printf(" fLongModuleSize %6.3f cm \n", fLongModuleSize);
424 printf(" #supermodule in phi direction %i \n", fNPhiSuperModule );
425 printf(" fILOSS %i : fIHADR %i \n", fILOSS, fIHADR);
426 printf(" fTrd1Angle %7.4f\n", fTrd1Angle);
427 printf(" f2Trd1Dx2 %7.4f\n", f2Trd1Dx2);
428 printf("SM dimensions(TRD1) : dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n",
429 fParSM[0],fParSM[1],fParSM[2]);
430 printf(" fPhiGapForSM %7.4f cm (%7.4f <- phi size in degree)\n",
431 fPhiGapForSM, TMath::ATan2(fPhiGapForSM,fIPDistance)*TMath::RadToDeg());
432 if(GetKey110DEG()) printf(" Last two modules have size 10 degree in phi (180<phi<190)\n");
433 printf(" phi SM boundaries \n");
434 for(int i=0; i<fPhiBoundariesOfSM.GetSize()/2.; i++) {
435 printf(" %i : %7.5f(%7.2f) -> %7.5f(%7.2f) : center %7.5f(%7.2f) \n", i,
436 fPhiBoundariesOfSM[2*i], fPhiBoundariesOfSM[2*i]*TMath::RadToDeg(),
437 fPhiBoundariesOfSM[2*i+1], fPhiBoundariesOfSM[2*i+1]*TMath::RadToDeg(),
438 fPhiCentersOfSM[i], fPhiCentersOfSM[i]*TMath::RadToDeg());
440 printf(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n",
441 fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1);
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");
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());
466 //______________________________________________________________________
467 void AliEMCALGeometry::PrintCellIndexes(Int_t absId, int pri, char *tit)
470 Int_t nSupMod, nModule, nIphi, nIeta;
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);
477 GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
478 printf(" local SM index : iphi %i : ieta %i \n", iphi,ieta);
479 GetGlobal(absId, vg);
480 printf(" vglob : mag %7.2f : perp %7.2f : z %7.2f : eta %6.4f : phi %6.4f(%6.2f) \n",
481 vg.Mag(), vg.Perp(), vg.Z(), vg.Eta(), vg.Phi(), vg.Phi()*TMath::RadToDeg());
485 //______________________________________________________________________
486 void AliEMCALGeometry::CheckAdditionalOptions()
489 // Additional options that
490 // can be used to select
491 // the specific geometry of
494 // adeed allILOSS= and allIHADR= for MIP investigation
495 fArrayOpts = new TObjArray;
496 Int_t nopt = AliEMCALHistoUtilities::ParseString(fGeoName, *fArrayOpts);
497 if(nopt==1) { // no aditional option(s)
498 fArrayOpts->Delete();
503 for(Int_t i=1; i<nopt; i++){
504 TObjString *o = (TObjString*)fArrayOpts->At(i);
506 TString addOpt = o->String();
508 for(Int_t j=0; j<fNAdditionalOpts; j++) {
509 TString opt = fAdditionalOpts[j];
510 if(addOpt.Contains(opt,TString::kIgnoreCase)) {
516 AliDebug(2,Form("<E> option |%s| unavailable : ** look to the file AliEMCALGeometry.h **\n",
520 AliDebug(2,Form("<I> option |%s| is valid : number %i : |%s|\n",
521 addOpt.Data(), indj, fAdditionalOpts[indj]));
522 if (addOpt.Contains("NL=",TString::kIgnoreCase)) {// number of sampling layers
523 sscanf(addOpt.Data(),"NL=%i", &fNECLayers);
524 AliDebug(2,Form(" fNECLayers %i (new) \n", fNECLayers));
525 } else if(addOpt.Contains("PBTH=",TString::kIgnoreCase)) {//Thickness of the Pb(fECPbRadThicknes)
526 sscanf(addOpt.Data(),"PBTH=%f", &fECPbRadThickness);
527 } else if(addOpt.Contains("SCTH=",TString::kIgnoreCase)) {//Thickness of the Sc(fECScintThick)
528 sscanf(addOpt.Data(),"SCTH=%f", &fECScintThick);
529 } else if(addOpt.Contains("LATSS=",TString::kIgnoreCase)) {// Thickness of lateral steel strip (fLateralSteelStrip)
530 sscanf(addOpt.Data(),"LATSS=%f", &fLateralSteelStrip);
531 AliDebug(2,Form(" fLateralSteelStrip %f (new) \n", fLateralSteelStrip));
532 } else if(addOpt.Contains("ILOSS=",TString::kIgnoreCase)) {// As in Geant
533 sscanf(addOpt.Data(),"ALLILOSS=%i", &fILOSS);
534 AliDebug(2,Form(" fILOSS %i \n", fILOSS));
535 } else if(addOpt.Contains("IHADR=",TString::kIgnoreCase)) {// As in Geant
536 sscanf(addOpt.Data(),"ALLIHADR=%i", &fIHADR);
537 AliDebug(2,Form(" fIHADR %i \n", fIHADR));
543 //__________________________________________________________________
544 void AliEMCALGeometry::DefineSamplingFraction()
547 // Look http://rhic.physics.wayne.edu/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
548 // Keep for compatibilty
550 if(fNECLayers == 69) { // 10% layer reduction
552 } else if(fNECLayers == 61) { // 20% layer reduction
554 } else if(fNECLayers == 77) {
555 if (fECScintThick>0.159 && fECScintThick<0.161) { // original sampling fraction, equal layers
556 fSampling = 12.327; // fECScintThick = fECPbRadThickness = 0.160;
557 } else if (fECScintThick>0.175 && fECScintThick<0.177) { // 10% Pb thicknes reduction
558 fSampling = 10.5; // fECScintThick = 0.176, fECPbRadThickness=0.144;
559 } else if(fECScintThick>0.191 && fECScintThick<0.193) { // 20% Pb thicknes reduction
560 fSampling = 8.93; // fECScintThick = 0.192, fECPbRadThickness=0.128;
566 //______________________________________________________________________
567 void AliEMCALGeometry::GetModulePhiEtaIndexInSModuleFromTRUIndex(const Int_t itru, const Int_t iphitru, const Int_t ietatru, Int_t &iphiSM, Int_t &ietaSM) const
570 // This method transforms the (eta,phi) index of module in a
571 // TRU matrix into Super Module (eta,phi) index.
573 // Calculate in which row and column where the TRU are
576 Int_t col = itru/ fNTRUPhi ; // indexes of TRU in SM
577 Int_t row = itru - col*fNTRUPhi ;
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);
585 //______________________________________________________________________
586 AliEMCALGeometry * AliEMCALGeometry::GetInstance(){
587 // Returns the pointer of the unique instance
589 AliEMCALGeometry * rv = static_cast<AliEMCALGeometry *>( fgGeom );
593 //______________________________________________________________________
594 AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
595 const Text_t* title){
596 // Returns the pointer of the unique instance
598 AliEMCALGeometry * rv = 0;
600 if ( strcmp(name,"") == 0 ) { // get default geometry
601 fgGeom = new AliEMCALGeometry(fgDefaultGeometryName, title);
603 fgGeom = new AliEMCALGeometry(name, title);
604 } // end if strcmp(name,"")
605 if ( fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
612 if ( strcmp(fgGeom->GetName(), name) != 0) {
613 printf("\ncurrent geometry is %s : ", fgGeom->GetName());
614 printf(" you cannot call %s ",name);
616 rv = (AliEMCALGeometry *) fgGeom;
622 //_____________________________________________________________________________
623 Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const {
624 // Checks whether point is inside the EMCal volume, used in AliEMCALv*.cxx
626 // Code uses cylindrical approximation made of inner radius (for speed)
628 // Points behind EMCAl, i.e. R > outer radius, but eta, phi in acceptance
629 // are considered to inside
631 Double_t r=sqrt(x*x+y*y);
633 if ( r > fEnvelop[0] ) {
635 theta = TMath::ATan2(r,z);
640 eta = -TMath::Log(TMath::Tan(theta/2.));
641 if (eta < fArm1EtaMin || eta > fArm1EtaMax)
644 Double_t phi = TMath::ATan2(y,x) * 180./TMath::Pi();
645 if (phi < 0) phi += 360; // phi should go from 0 to 360 in this case
646 if (phi > fArm1PhiMin && phi < fArm1PhiMax)
653 // == Shish-kebab cases ==
655 //________________________________________________________________________________________________
656 Int_t AliEMCALGeometry::GetAbsCellId(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta) const
660 // 13-oct-05; 110 degree case
661 // May 31, 2006; ALICE numbering scheme:
662 // 0 <= nSupMod < fNumberOfSuperModules
663 // 0 <= nModule < fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1)
664 // 0 <= nIphi < fNPHIdiv
665 // 0 <= nIeta < fNETAdiv
666 // 0 <= absid < fNCells
667 static Int_t id=0; // have to change from 0 to fNCells-1
668 if(fKey110DEG == 1 && nSupMod >= 10) { // 110 degree case; last two supermodules
669 id = fNCellsInSupMod*10 + (fNCellsInSupMod/2)*(nSupMod-10);
671 id = fNCellsInSupMod*nSupMod;
673 id += fNCellsInModule *nModule;
674 id += fNPHIdiv *nIphi;
676 if(id<0 || id >= fNCells) {
677 // printf(" wrong numerations !!\n");
678 // printf(" id %6i(will be force to -1)\n", id);
679 // printf(" fNCells %6i\n", fNCells);
680 // printf(" nSupMod %6i\n", nSupMod);
681 // printf(" nModule %6i\n", nModule);
682 // printf(" nIphi %6i\n", nIphi);
683 // printf(" nIeta %6i\n", nIeta);
684 id = -TMath::Abs(id); // if negative something wrong
689 //________________________________________________________________________________________________
690 Bool_t AliEMCALGeometry::CheckAbsCellId(Int_t absId) const
692 // May 31, 2006; only trd1 now
693 if(absId<0 || absId >= fNCells) return kFALSE;
697 //________________________________________________________________________________________________
698 Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nModule,Int_t &nIphi,Int_t &nIeta) const
700 // 21-sep-04; 19-oct-05;
701 // May 31, 2006; ALICE numbering scheme:
704 // absId - cell is as in Geant, 0<= absId < fNCells;
706 // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
707 // nModule - module number in SM, 0<= nModule < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
708 // nIphi - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv;
709 // nIeta - cell number in eta driection inside module; 0<= nIeta < fNETAdiv;
711 static Int_t tmp=0, sm10=0;
712 if(!CheckAbsCellId(absId)) return kFALSE;
714 sm10 = fNCellsInSupMod*10;
715 if(fKey110DEG == 1 && absId >= sm10) { // 110 degree case; last two supermodules
716 nSupMod = (absId-sm10) / (fNCellsInSupMod/2) + 10;
717 tmp = (absId-sm10) % (fNCellsInSupMod/2);
719 nSupMod = absId / fNCellsInSupMod;
720 tmp = absId % fNCellsInSupMod;
723 nModule = tmp / fNCellsInModule;
724 tmp = tmp % fNCellsInModule;
725 nIphi = tmp / fNPHIdiv;
726 nIeta = tmp % fNPHIdiv;
731 //________________________________________________________________________________________________
732 void AliEMCALGeometry::GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, int &iphim, int &ietam) const
734 // added nSupMod; - 19-oct-05 !
735 // Alice numbering scheme - Jun 01,2006
736 // ietam, iphi - indexes of module in two dimensional grid of SM
737 // ietam - have to change from 0 to fNZ-1
738 // iphim - have to change from 0 to nphi-1 (fNPhi-1 or fNPhi/2-1)
741 if(fKey110DEG == 1 && nSupMod>=10) nphi = fNPhi/2;
744 ietam = nModule/nphi;
745 iphim = nModule%nphi;
748 //________________________________________________________________________________________________
749 void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta,
750 int &iphi, int &ieta) const
753 // Added nSupMod; Nov 25, 05
754 // Alice numbering scheme - Jun 01,2006
756 // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
757 // nModule - module number in SM, 0<= nModule < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
758 // nIphi - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv;
759 // nIeta - cell number in eta driection inside module; 0<= nIeta < fNETAdiv;
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)
766 static Int_t iphim, ietam;
768 GetModulePhiEtaIndexInSModule(nSupMod,nModule, iphim, ietam);
769 // ieta = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM)
770 ieta = ietam*fNETAdiv + (fNETAdiv - 1 - nIeta); // x(module) = -z(SM)
771 iphi = iphim*fNPHIdiv + nIphi; // y(module) = y(SM)
774 AliDebug(1,Form(" nSupMod %i nModule %i nIphi %i nIeta %i => ieta %i iphi %i\n",
775 nSupMod, nModule, nIphi, nIeta, ieta, iphi));
778 //________________________________________________________________________________________________
779 Int_t AliEMCALGeometry::GetSuperModuleNumber(Int_t absId) const
781 // Return the number of the supermodule given the absolute
782 // ALICE numbering id
784 static Int_t nSupMod, nModule, nIphi, nIeta;
785 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
789 //________________________________________________________________________________________________
790 void AliEMCALGeometry::GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta,
791 Int_t &iphim, Int_t &ietam, Int_t &nModule) const
793 // Transition from cell indexes (ieta,iphi) to module indexes (ietam,iphim, nModule)
795 nphi = GetNumberOfModuleInPhiDirection(nSupMod);
797 ietam = ieta/fNETAdiv;
798 iphim = iphi/fNPHIdiv;
799 nModule = ietam * nphi + iphim;
802 //________________________________________________________________________________________________
803 Int_t AliEMCALGeometry::GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const
805 // Transition from super module number(nSupMod) and cell indexes (ieta,iphi) to absId
806 static Int_t ietam, iphim, nModule;
807 static Int_t nIeta, nIphi; // cell indexes in module
809 GetModuleIndexesFromCellIndexesInSModule(nSupMod, iphi, ieta, ietam, iphim, nModule);
811 nIeta = ieta%fNETAdiv;
812 nIeta = fNETAdiv - 1 - nIeta;
813 nIphi = iphi%fNPHIdiv;
815 return GetAbsCellId(nSupMod, nModule, nIphi, nIeta);
819 // Methods for AliEMCALRecPoint - Feb 19, 2006
820 //________________________________________________________________________________________________
821 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const
823 // Look to see what the relative
824 // position inside a given cell is
826 // Alice numbering scheme - Jun 08, 2006
828 // absId - cell is as in Geant, 0<= absId < fNCells;
830 // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
832 // Shift index taking into account the difference between standard SM
833 // and SM of half size in phi direction
834 const Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
835 static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
836 if(!CheckAbsCellId(absId)) return kFALSE;
838 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
839 GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
841 xr = fCentersOfCellsXDir.At(ieta);
842 zr = fCentersOfCellsEtaDir.At(ieta);
845 yr = fCentersOfCellsPhiDir.At(iphi);
847 yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
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));
854 //________________________________________________________________________________________________
855 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t loc[3]) const
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])) {
865 //________________________________________________________________________________________________
866 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, TVector3 &vloc) const
868 static Double_t loc[3];
869 if(RelPosCellInSModule(absId,loc)) {
870 vloc.SetXYZ(loc[0], loc[1], loc[2]);
876 // Alice numbering scheme - Jun 03, 2006
879 //________________________________________________________________________________________________
880 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t distEff, Double_t &xr, Double_t &yr, Double_t &zr) const
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
887 // absId - cell is as in Geant, 0<= absId < fNCells;
888 // e - cluster energy
890 // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
892 // Shift index taking into account the difference between standard SM
893 // and SM of half size in phi direction
894 const Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
895 static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
896 static Int_t iphim, ietam;
897 static AliEMCALShishKebabTrd1Module *mod = 0;
899 if(!CheckAbsCellId(absId)) return kFALSE;
901 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
902 GetModulePhiEtaIndexInSModule(nSupMod, nModule, iphim, ietam);
903 GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
905 mod = GetShishKebabModule(ietam);
906 mod->GetPositionAtCenterCellLine(nIeta, distEff, v);
907 xr = v.Y() - fParSM[0];
908 zr = v.X() - fParSM[2];
911 yr = fCentersOfCellsPhiDir.At(iphi);
913 yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
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));
920 //________________________________________________________________________________________________
921 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Int_t maxAbsId, Double_t distEff, Double_t &xr, Double_t &yr, Double_t &zr) const
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
928 // absId - cell is as in Geant, 0<= absId < fNCells;
929 // maxAbsId - abs id of cell with highest energy
930 // e - cluster energy
932 // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
934 // Shift index taking into account the difference between standard SM
935 // and SM of half size in phi direction
936 const Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
937 static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
938 static Int_t iphim, ietam;
939 static AliEMCALShishKebabTrd1Module *mod = 0;
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;
947 if(!CheckAbsCellId(absId)) return kFALSE;
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);
954 if(absId != maxAbsId) {
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;
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);
967 // distEff += distCorr;
969 // Bad resolution in this case, strong bias vs phi
971 mod->GetPositionAtCenterCellLine(nIeta, distEff, v); // Stay here
972 xr = v.Y() - fParSM[0];
973 zr = v.X() - fParSM[2];
976 yr = fCentersOfCellsPhiDir.At(iphi);
978 yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
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));
985 //________________________________________________________________________________________________
986 void AliEMCALGeometry::CreateListOfTrd1Modules()
988 // Generate the list of Trd1 modules
989 // which will make up the EMCAL
992 AliDebug(2,Form(" AliEMCALGeometry::CreateListOfTrd1Modules() started "));
994 AliEMCALShishKebabTrd1Module *mod=0, *mTmp=0; // current module
995 if(fShishKebabTrd1Modules == 0) {
996 fShishKebabTrd1Modules = new TList;
997 fShishKebabTrd1Modules->SetName("ListOfTRD1");
998 for(int iz=0; iz< GetNZ(); iz++) {
1000 mod = new AliEMCALShishKebabTrd1Module(TMath::Pi()/2.,this);
1002 mTmp = new AliEMCALShishKebabTrd1Module(*mod);
1005 fShishKebabTrd1Modules->Add(mod);
1008 AliDebug(2,Form(" Already exits : "));
1010 mod = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(fShishKebabTrd1Modules->GetSize()-1);
1011 fEtaMaxOfTRD1 = mod->GetMaxEtaOfModule(0);
1013 AliDebug(2,Form(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n",
1014 fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1));
1016 // Jun 01, 2006 - ALICE numbering scheme
1017 // define grid for cells in eta(z) and x directions in local coordinates system of SM
1018 // Works just for 2x2 case only -- ?? start here
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
1024 AliDebug(2,Form(" Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize()));
1025 Int_t ind=0; // this is phi index
1026 Int_t ieta=0, nModule=0, iphiTemp;
1027 Double_t xr, zr, theta, phi, eta, r, x,y;
1029 Double_t ytCenterModule=0.0, ytCenterCell=0.0;
1031 fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
1032 fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
1034 Double_t r0 = GetIPDistance() + GetLongModuleSize()/2.;
1035 for(Int_t it=0; it<fNPhi; it++) { // cycle on modules
1036 ytCenterModule = -fParSM[1] + fPhiModuleSize*(2*it+1)/2; // center of module
1037 for(Int_t ic=0; ic<fNPHIdiv; ic++) { // cycle on cells in module
1039 ytCenterCell = ytCenterModule + fPhiTileSize *(2*ic-1)/2.;
1040 } else if(fNPHIdiv==3){
1041 ytCenterCell = ytCenterModule + fPhiTileSize *(ic-1);
1042 } else if(fNPHIdiv==1){
1043 ytCenterCell = ytCenterModule;
1045 fCentersOfCellsPhiDir.AddAt(ytCenterCell,ind);
1046 // Define grid on phi direction
1047 // Grid is not the same for different eta bin;
1048 // Effect is small but is still here
1049 phi = TMath::ATan2(ytCenterCell, r0);
1050 fPhiCentersOfCells.AddAt(phi, ind);
1052 AliDebug(2,Form(" ind %2.2i : y %8.3f ", ind, fCentersOfCellsPhiDir.At(ind)));
1057 fCentersOfCellsEtaDir.Set(fNZ *fNETAdiv);
1058 fCentersOfCellsXDir.Set(fNZ *fNETAdiv);
1059 fEtaCentersOfCells.Set(fNZ *fNETAdiv * fNPhi*fNPHIdiv);
1060 AliDebug(2,Form(" Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize()));
1061 for(Int_t it=0; it<fNZ; it++) {
1062 AliEMCALShishKebabTrd1Module *trd1 = GetShishKebabModule(it);
1064 for(Int_t ic=0; ic<fNETAdiv; ic++) {
1066 trd1->GetCenterOfCellInLocalCoordinateofSM(ic, xr, zr); // case of 2X2
1067 GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1069 trd1->GetCenterOfCellInLocalCoordinateofSM_3X3(ic, xr, zr); // case of 3X3
1070 GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1072 trd1->GetCenterOfCellInLocalCoordinateofSM_1X1(xr, zr); // case of 1X1
1073 GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1075 fCentersOfCellsXDir.AddAt(float(xr) - fParSM[0],ieta);
1076 fCentersOfCellsEtaDir.AddAt(float(zr) - fParSM[2],ieta);
1077 // Define grid on eta direction for each bin in phi
1078 for(int iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
1079 x = xr + trd1->GetRadius();
1080 y = fCentersOfCellsPhiDir[iphi];
1081 r = TMath::Sqrt(x*x + y*y + zr*zr);
1082 theta = TMath::ACos(zr/r);
1083 eta = AliEMCALShishKebabTrd1Module::ThetaToEta(theta);
1084 // ind = ieta*fCentersOfCellsPhiDir.GetSize() + iphi;
1085 ind = iphi*fCentersOfCellsEtaDir.GetSize() + ieta;
1086 fEtaCentersOfCells.AddAt(eta, ind);
1088 //printf(" ieta %i : xr + trd1->GetRadius() %f : zr %f : eta %f \n", ieta, xr + trd1->GetRadius(), zr, eta);
1091 for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1092 AliDebug(2,Form(" ind %2.2i : z %8.3f : x %8.3f", i+1,
1093 fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i)));
1098 //________________________________________________________________________________________________
1099 void AliEMCALGeometry::GetGlobal(const Double_t *loc, Double_t *glob, int ind) const
1101 // Figure out the global numbering
1102 // of a given supermodule from the
1103 // local numbering and the transformation
1104 // matrix stored by the geometry manager (allows for misaligned
1107 if(ind>=0 && ind < GetNumberOfSuperModules()) {
1108 TString volpath = "ALIC_1/XEN1_1/SMOD_";
1111 if(GetKey110DEG() && ind>=10) {
1112 volpath = "ALIC_1/XEN1_1/SM10_";
1113 volpath += ind-10+1;
1116 if(!gGeoManager->cd(volpath.Data()))
1117 AliFatal(Form("AliEMCALGeometry::GeoManager cannot find path %s!",volpath.Data()));
1119 TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1121 m->LocalToMaster(loc, glob);
1123 AliFatal("Geo matrixes are not loaded \n") ;
1128 //________________________________________________________________________________________________
1129 void AliEMCALGeometry::GetGlobal(const TVector3 &vloc, TVector3 &vglob, int ind) const
1131 //Figure out the global numbering
1132 //of a given supermodule from the
1133 //local numbering given a 3-vector location
1135 static Double_t tglob[3], tloc[3];
1137 GetGlobal(tloc, tglob, ind);
1138 vglob.SetXYZ(tglob[0], tglob[1], tglob[2]);
1141 //________________________________________________________________________________________________
1142 void AliEMCALGeometry::GetGlobal(Int_t absId , double glob[3]) const
1144 // Alice numbering scheme - Jun 03, 2006
1145 static Int_t nSupMod, nModule, nIphi, nIeta;
1146 static double loc[3];
1148 if (!gGeoManager || !gGeoManager->IsClosed()) {
1149 AliError("Can't get the global coordinates! gGeoManager doesn't exist or it is still open!");
1153 glob[0]=glob[1]=glob[2]=0.0; // bad case
1154 if(RelPosCellInSModule(absId, loc)) {
1155 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1157 TString volpath = "ALIC_1/XEN1_1/SMOD_";
1158 volpath += (nSupMod+1);
1160 if(GetKey110DEG() && nSupMod>=10) {
1161 volpath = "ALIC_1/XEN1_1/SM10_";
1162 volpath += (nSupMod-10+1);
1164 if(!gGeoManager->cd(volpath.Data()))
1165 AliFatal(Form("GeoManager cannot find path %s!",volpath.Data()));
1167 TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1169 m->LocalToMaster(loc, glob);
1171 AliFatal("Geo matrixes are not loaded \n") ;
1176 //___________________________________________________________________
1177 void AliEMCALGeometry::GetGlobal(Int_t absId , TVector3 &vglob) const
1179 // Alice numbering scheme - Jun 03, 2006
1180 static Double_t glob[3];
1182 GetGlobal(absId, glob);
1183 vglob.SetXYZ(glob[0], glob[1], glob[2]);
1187 //____________________________________________________________________________
1188 void AliEMCALGeometry::GetGlobal(const AliRecPoint* /*rp*/, TVector3& /* vglob */) const
1190 AliFatal(Form("Please use GetGlobalEMCAL(recPoint,gpos) instead of GetGlobal!"));
1193 //_________________________________________________________________________________
1194 void AliEMCALGeometry::GetGlobalEMCAL(const AliEMCALRecPoint *rp, TVector3 &vglob) const
1196 // Figure out the global numbering
1197 // of a given supermodule from the
1198 // local numbering for RecPoints
1200 static TVector3 vloc;
1201 static Int_t nSupMod, nModule, nIphi, nIeta;
1203 const AliEMCALRecPoint *rpTmp = rp;
1204 const AliEMCALRecPoint *rpEmc = rpTmp;
1206 GetCellIndex(rpEmc->GetAbsId(0), nSupMod, nModule, nIphi, nIeta);
1207 rpTmp->GetLocalPosition(vloc);
1208 GetGlobal(vloc, vglob, nSupMod);
1211 //________________________________________________________________________________________________
1212 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Double_t &eta,Double_t &phi) const
1214 // Nov 16, 2006- float to double
1215 // version for TRD1 only
1216 static TVector3 vglob;
1217 GetGlobal(absId, vglob);
1222 //________________________________________________________________________________________________
1223 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Float_t &eta,Float_t &phi) const
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());
1232 //________________________________________________________________________________________________
1233 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSM(Int_t nSupMod, Double_t &phiMin, Double_t &phiMax) const
1235 // 0<= nSupMod <=11; phi in rad
1237 if(nSupMod<0 || nSupMod >11) return kFALSE;
1239 phiMin = fPhiBoundariesOfSM[2*i];
1240 phiMax = fPhiBoundariesOfSM[2*i+1];
1244 //________________________________________________________________________________________________
1245 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const
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];
1259 //________________________________________________________________________________________________
1260 Bool_t AliEMCALGeometry::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const
1262 // Return false if phi belongs a phi cracks between SM
1266 if(TMath::Abs(eta) > fEtaMaxOfTRD1) return kFALSE;
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]) {
1272 if(eta < 0.0) nSupMod++;
1273 AliDebug(1,Form("eta %f phi %f(%5.2f) : nSupMod %i : #bound %i", eta,phi,phi*TMath::RadToDeg(), nSupMod,i));
1280 //________________________________________________________________________________________________
1281 Bool_t AliEMCALGeometry::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi, Int_t &absId) const
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)) {
1290 phi = TVector2::Phi_0_2pi(phi);
1291 phiLoc = phi - fPhiCentersOfSM[nSupMod/2];
1292 nphi = fPhiCentersOfCells.GetSize();
1294 phiLoc = phi - 190.*TMath::DegToRad();
1298 dmin = TMath::Abs(fPhiCentersOfCells[0]-phiLoc);
1300 for(i=1; i<nphi; i++) {
1301 d = TMath::Abs(fPhiCentersOfCells[i] - phiLoc);
1306 // printf(" i %i : d %f : dmin %f : fPhiCentersOfCells[i] %f \n", i, d, dmin, fPhiCentersOfCells[i]);
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));
1311 absEta = TMath::Abs(eta);
1312 etaShift = iphi*fCentersOfCellsEtaDir.GetSize();
1313 dmin = TMath::Abs(fEtaCentersOfCells[etaShift]-absEta);
1315 for(i=1; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1316 d = TMath::Abs(fEtaCentersOfCells[i+etaShift] - absEta);
1322 AliDebug(2,Form(" ieta %i : dmin %f (eta=%f) : nSupMod %i ", ieta, dmin, eta, nSupMod));
1324 if(eta<0) iphi = (nphi-1) - iphi;
1325 absId = GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);
1332 //________________________________________________________________________________________________
1333 AliEMCALShishKebabTrd1Module* AliEMCALGeometry::GetShishKebabModule(Int_t neta) const
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.
1341 static AliEMCALShishKebabTrd1Module* trd1=0;
1342 if(fShishKebabTrd1Modules && neta>=0 && neta<fShishKebabTrd1Modules->GetSize()) {
1343 trd1 = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(neta);
1348 //________________________________________________________________________________________________
1349 Int_t AliEMCALGeometry::GetAbsTRUNumberFromNumberInSm(const Int_t row, const Int_t col, const Int_t sm)
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);
1356 //________________________________________________________________________________________________
1357 void AliEMCALGeometry::Browse(TBrowser* b)
1359 //Browse the modules
1360 if(fShishKebabTrd1Modules) b->Add(fShishKebabTrd1Modules);
1363 //________________________________________________________________________________________________
1364 Bool_t AliEMCALGeometry::IsFolder() const
1366 //Check if fShishKebabTrd1Modules is in folder
1367 if(fShishKebabTrd1Modules) return kTRUE;
1371 //________________________________________________________________________________________________
1372 Double_t AliEMCALGeometry::GetPhiCenterOfSM(Int_t nsupmod) const
1374 //returns center of supermodule in phi
1376 return fPhiCentersOfSM[i];
1379 //____________________________________________________________________________
1380 Bool_t AliEMCALGeometry::Impact(const TParticle * particle) const
1382 // Tells if a particle enters EMCAL
1385 TVector3 vtx(particle->Vx(),particle->Vy(),particle->Vz());
1386 TVector3 vimpact(0,0,0);
1387 ImpactOnEmcal(vtx,particle->Theta(),particle->Phi(),AbsID,vimpact);
1392 //____________________________________________________________________________
1393 void AliEMCALGeometry::ImpactOnEmcal(TVector3 vtx, Double_t theta, Double_t phi,
1394 Int_t & absId, TVector3 & vimpact) const
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
1400 TVector3 p(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) ;
1402 vimpact.SetXYZ(0,0,0);
1404 if(phi==0 || theta==0) return;
1407 Double_t factor = (GetIPDistance()-vtx[1])/p[1];
1408 direction = vtx + factor*p;
1411 AliFatal("Geo manager not initialized\n");
1413 //from particle direction -> tower hitted
1414 GetAbsCellIdFromEtaPhi(direction.Eta(),direction.Phi(),absId);
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]={};
1421 if(!RelPosCellInSModule(absId,loc)) return;
1423 //loc is cell center of tower
1424 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
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;
1430 absId2=GetAbsCellId(nSupMod,nModule,nIphi,nIeta2);
1431 if(nIphi==0) nIphi2=1;
1433 absId3=GetAbsCellId(nSupMod,nModule,nIphi2,nIeta);
1435 //2nd point on emcal cell plane
1436 if(!RelPosCellInSModule(absId2,loc2)) return;
1438 //3rd point on emcal cell plane
1439 if(!RelPosCellInSModule(absId3,loc3)) return;
1441 TString volpath = "ALIC_1/XEN1_1/SMOD_";
1442 volpath += (nSupMod+1);
1444 if(GetKey110DEG() && nSupMod>=10) {
1445 volpath = "ALIC_1/XEN1_1/SM10_";
1446 volpath += (nSupMod-10+1);
1448 if(!gGeoManager->cd(volpath.Data())){
1449 AliFatal(Form("GeoManager cannot find path %s!",volpath.Data()))
1452 TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1454 m->LocalToMaster(loc, glob);
1455 m->LocalToMaster(loc2, glob2);
1456 m->LocalToMaster(loc3, glob3);
1458 AliFatal("Geo matrixes are not loaded \n") ;
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]);
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 ;
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));
1485 printf("ImpactOnEmcal() No solution :\n");
1489 Double_t length = A*vtx(0)+B*vtx(1)+C*vtx(2)+D;
1492 vimpact.SetXYZ(vtx(0)+length*(direction(0)-vtx(0)),vtx(1)+length*(direction(1)-vtx(1)),vtx(2)+length*(direction(2)-vtx(2)));
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);