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 const 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::Instance();
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") || fGeoName.Contains("EMCAL_1stYear"))) {
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 if(fGeoName.Contains("1stYear")){
335 fNumberOfSuperModules = 2;
337 if(fGeoName.Contains("LowerEta")) {
338 fNPhiSuperModule = 1;
340 else if(fGeoName.Contains("LowerPhi_SideA")){
341 fNPhiSuperModule = 2;
344 else if(fGeoName.Contains("LowerPhi_SideC")){
345 fNPhiSuperModule = 2;
349 CheckAdditionalOptions();
352 // constant for transition absid <--> indexes
353 fNCellsInModule = fNPHIdiv*fNETAdiv;
354 fNCellsInSupMod = fNCellsInModule*fNPhi*fNZ;
355 fNCells = fNCellsInSupMod*fNumberOfSuperModules;
356 if(GetKey110DEG()) fNCells -= fNCellsInSupMod;
358 fNPhiSuperModule = fNumberOfSuperModules/2;
359 if(fNPhiSuperModule < 1) fNPhiSuperModule = 1;
361 fPhiTileSize = fPhiModuleSize/double(fNPHIdiv) - fLateralSteelStrip; // 13-may-05
362 fEtaTileSize = fEtaModuleSize/double(fNETAdiv) - fLateralSteelStrip; // 13-may-05
364 fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness);
365 f2Trd1Dx2 = fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.);
366 if(!fGeoName.Contains("WSUC")) fShellThickness = TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2);
368 //These parameters are used to create the mother volume to hold the supermodules
369 //2cm padding added to allow for misalignments - JLK 30-May-2008
370 fEnvelop[0] = fIPDistance - 1.; // mother volume inner radius
371 fEnvelop[1] = fIPDistance + fShellThickness + 1.; // mother volume outer r.
372 fEnvelop[2] = fZLength + 2.; //mother volume length
375 fParSM[0] = GetShellThickness()/2.;
376 fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
377 fParSM[2] = fZLength/4.; //divide by 4 to get half-length of SM
379 // SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006
380 fPhiBoundariesOfSM.Set(fNumberOfSuperModules);
381 fPhiCentersOfSM.Set(fNumberOfSuperModules/2);
382 fPhiBoundariesOfSM[0] = TMath::PiOver2() - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules)
383 fPhiCentersOfSM[0] = TMath::PiOver2();
384 if(fNumberOfSuperModules > 1)
385 fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance);
386 if(fNumberOfSuperModules > 2) {
387 for(int i=1; i<=4; i++) { // from 2th ro 9th
388 fPhiBoundariesOfSM[2*i] = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i;
389 fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i;
390 fPhiCentersOfSM[i] = fPhiCentersOfSM[0] + 20.*TMath::DegToRad()*i;
393 if(fNumberOfSuperModules > 10) {
394 fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad();
395 fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance);
396 fPhiCentersOfSM[5] = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.;
399 //called after setting of scintillator and lead layer parameters
400 DefineSamplingFraction();
402 // TRU parameters - Apr 29,08 by PAI.
403 // These parameters values was updated at Nov 05, 2007
404 // As is on Olivier BOURRION (LPSC) ppt preasentation
405 // at ALICE trigger meeting at 13th-14th March
406 fNTRUEta = 1; // was 3
407 fNTRUPhi = 3; // was 1
408 fNModulesInTRUEta = 24; // was 8
409 fNModulesInTRUPhi = 4; // was 12
411 // 3*6*10 + 2*6*2 = 204 -> matrix (nphi(17), neta(12))
417 //___________________________________________________________________
418 void AliEMCALGeometry::PrintGeometry()
420 // Separate routine is callable from broswer; Nov 7,2006
421 printf("\nInit: geometry of EMCAL named %s :\n", fGeoName.Data());
423 for(Int_t i=0; i<fArrayOpts->GetEntries(); i++){
424 TObjString *o = (TObjString*)fArrayOpts->At(i);
425 printf(" %i : %s \n", i, o->String().Data());
428 printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
429 printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f -> for EMCAL envelope only\n",
430 GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() );
432 printf( " ECAL : %d x (%f cm Pb, %f cm Sc) \n",
433 GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ;
434 printf(" fSampling %5.2f \n", fSampling );
435 printf(" fIPDistance %6.3f cm \n", fIPDistance);
436 printf(" fNPhi %i | fNZ %i \n", fNPhi, fNZ);
437 printf(" fNCellsInModule %i : fNCellsInSupMod %i : fNCells %i\n",fNCellsInModule, fNCellsInSupMod, fNCells);
438 printf(" X:Y module size %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize);
439 printf(" X:Y tile size %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize);
440 printf(" #of sampling layers %i(fNECLayers) \n", fNECLayers);
441 printf(" fLongModuleSize %6.3f cm \n", fLongModuleSize);
442 printf(" #supermodule in phi direction %i \n", fNPhiSuperModule );
443 printf(" fILOSS %i : fIHADR %i \n", fILOSS, fIHADR);
444 printf(" fTrd1Angle %7.4f\n", fTrd1Angle);
445 printf(" f2Trd1Dx2 %7.4f\n", f2Trd1Dx2);
446 printf("SM dimensions(TRD1) : dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n",
447 fParSM[0],fParSM[1],fParSM[2]);
448 printf(" fPhiGapForSM %7.4f cm (%7.4f <- phi size in degree)\n",
449 fPhiGapForSM, TMath::ATan2(fPhiGapForSM,fIPDistance)*TMath::RadToDeg());
450 if(GetKey110DEG()) printf(" Last two modules have size 10 degree in phi (180<phi<190)\n");
451 printf(" phi SM boundaries \n");
452 for(int i=0; i<fPhiBoundariesOfSM.GetSize()/2.; i++) {
453 printf(" %i : %7.5f(%7.2f) -> %7.5f(%7.2f) : center %7.5f(%7.2f) \n", i,
454 fPhiBoundariesOfSM[2*i], fPhiBoundariesOfSM[2*i]*TMath::RadToDeg(),
455 fPhiBoundariesOfSM[2*i+1], fPhiBoundariesOfSM[2*i+1]*TMath::RadToDeg(),
456 fPhiCentersOfSM[i], fPhiCentersOfSM[i]*TMath::RadToDeg());
458 printf(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n",
459 fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1);
461 printf("\n Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize());
462 for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
463 printf(" ind %2.2i : z %8.3f : x %8.3f \n", i,
464 fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i));
465 int ind=0; // Nov 21,2006
466 for(Int_t iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
467 ind = iphi*fCentersOfCellsEtaDir.GetSize() + i;
468 printf("%6.4f ", fEtaCentersOfCells[ind]);
469 if((iphi+1)%12 == 0) printf("\n");
475 printf("\n Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize());
476 for(Int_t i=0; i<fCentersOfCellsPhiDir.GetSize(); i++) {
477 double phi=fPhiCentersOfCells.At(i);
478 printf(" ind %2.2i : y %8.3f : phi %7.5f(%6.2f) \n", i, fCentersOfCellsPhiDir.At(i),
479 phi, phi*TMath::RadToDeg());
484 //______________________________________________________________________
485 void AliEMCALGeometry::PrintCellIndexes(Int_t absId, int pri, const char *tit)
488 Int_t nSupMod, nModule, nIphi, nIeta;
492 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
493 printf(" %s | absId : %i -> nSupMod %i nModule %i nIphi %i nIeta %i \n", tit, absId, nSupMod, nModule, nIphi, nIeta);
495 GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
496 printf(" local SM index : iphi %i : ieta %i \n", iphi,ieta);
497 GetGlobal(absId, vg);
498 printf(" vglob : mag %7.2f : perp %7.2f : z %7.2f : eta %6.4f : phi %6.4f(%6.2f) \n",
499 vg.Mag(), vg.Perp(), vg.Z(), vg.Eta(), vg.Phi(), vg.Phi()*TMath::RadToDeg());
503 //______________________________________________________________________
504 void AliEMCALGeometry::CheckAdditionalOptions()
507 // Additional options that
508 // can be used to select
509 // the specific geometry of
512 // adeed allILOSS= and allIHADR= for MIP investigation
513 fArrayOpts = new TObjArray;
514 Int_t nopt = AliEMCALHistoUtilities::ParseString(fGeoName, *fArrayOpts);
515 if(nopt==1) { // no aditional option(s)
516 fArrayOpts->Delete();
521 for(Int_t i=1; i<nopt; i++){
522 TObjString *o = (TObjString*)fArrayOpts->At(i);
524 TString addOpt = o->String();
526 for(Int_t j=0; j<fNAdditionalOpts; j++) {
527 TString opt = fAdditionalOpts[j];
528 if(addOpt.Contains(opt,TString::kIgnoreCase)) {
534 AliDebug(2,Form("<E> option |%s| unavailable : ** look to the file AliEMCALGeometry.h **\n",
538 AliDebug(2,Form("<I> option |%s| is valid : number %i : |%s|\n",
539 addOpt.Data(), indj, fAdditionalOpts[indj]));
540 if (addOpt.Contains("NL=",TString::kIgnoreCase)) {// number of sampling layers
541 sscanf(addOpt.Data(),"NL=%i", &fNECLayers);
542 AliDebug(2,Form(" fNECLayers %i (new) \n", fNECLayers));
543 } else if(addOpt.Contains("PBTH=",TString::kIgnoreCase)) {//Thickness of the Pb(fECPbRadThicknes)
544 sscanf(addOpt.Data(),"PBTH=%f", &fECPbRadThickness);
545 } else if(addOpt.Contains("SCTH=",TString::kIgnoreCase)) {//Thickness of the Sc(fECScintThick)
546 sscanf(addOpt.Data(),"SCTH=%f", &fECScintThick);
547 } else if(addOpt.Contains("LATSS=",TString::kIgnoreCase)) {// Thickness of lateral steel strip (fLateralSteelStrip)
548 sscanf(addOpt.Data(),"LATSS=%f", &fLateralSteelStrip);
549 AliDebug(2,Form(" fLateralSteelStrip %f (new) \n", fLateralSteelStrip));
550 } else if(addOpt.Contains("ILOSS=",TString::kIgnoreCase)) {// As in Geant
551 sscanf(addOpt.Data(),"ALLILOSS=%i", &fILOSS);
552 AliDebug(2,Form(" fILOSS %i \n", fILOSS));
553 } else if(addOpt.Contains("IHADR=",TString::kIgnoreCase)) {// As in Geant
554 sscanf(addOpt.Data(),"ALLIHADR=%i", &fIHADR);
555 AliDebug(2,Form(" fIHADR %i \n", fIHADR));
561 //__________________________________________________________________
562 void AliEMCALGeometry::DefineSamplingFraction()
565 // Look http://rhic.physics.wayne.edu/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
566 // Keep for compatibilty
568 if(fNECLayers == 69) { // 10% layer reduction
570 } else if(fNECLayers == 61) { // 20% layer reduction
572 } else if(fNECLayers == 77) {
573 if (fECScintThick>0.159 && fECScintThick<0.161) { // original sampling fraction, equal layers
574 fSampling = 12.327; // fECScintThick = fECPbRadThickness = 0.160;
575 } else if (fECScintThick>0.175 && fECScintThick<0.177) { // 10% Pb thicknes reduction
576 fSampling = 10.5; // fECScintThick = 0.176, fECPbRadThickness=0.144;
577 } else if(fECScintThick>0.191 && fECScintThick<0.193) { // 20% Pb thicknes reduction
578 fSampling = 8.93; // fECScintThick = 0.192, fECPbRadThickness=0.128;
584 //______________________________________________________________________
585 void AliEMCALGeometry::GetModulePhiEtaIndexInSModuleFromTRUIndex(Int_t itru, Int_t iphitru, Int_t ietatru, Int_t &iphiSM, Int_t &ietaSM) const
588 // This method transforms the (eta,phi) index of module in a
589 // TRU matrix into Super Module (eta,phi) index.
591 // Calculate in which row and column where the TRU are
594 Int_t col = itru/ fNTRUPhi ; // indexes of TRU in SM
595 Int_t row = itru - col*fNTRUPhi ;
597 iphiSM = fNModulesInTRUPhi*row + iphitru ;
598 ietaSM = fNModulesInTRUEta*col + ietatru ;
599 //printf(" GetModulePhiEtaIndexInSModuleFromTRUIndex : itru %2i iphitru %2i ietatru %2i iphiSM %2i ietaSM %2i \n",
600 // itru, iphitru, ietatru, iphiSM, ietaSM);
603 //______________________________________________________________________
604 AliEMCALGeometry * AliEMCALGeometry::GetInstance(){
605 // Returns the pointer of the unique instance
607 AliEMCALGeometry * rv = static_cast<AliEMCALGeometry *>( fgGeom );
611 //______________________________________________________________________
612 AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
613 const Text_t* title){
614 // Returns the pointer of the unique instance
616 AliEMCALGeometry * rv = 0;
618 if ( strcmp(name,"") == 0 ) { // get default geometry
619 fgGeom = new AliEMCALGeometry(fgDefaultGeometryName, title);
621 fgGeom = new AliEMCALGeometry(name, title);
622 } // end if strcmp(name,"")
623 if ( fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
630 if ( strcmp(fgGeom->GetName(), name) != 0) {
631 printf("\ncurrent geometry is %s : ", fgGeom->GetName());
632 printf(" you cannot call %s ",name);
634 rv = (AliEMCALGeometry *) fgGeom;
640 //_____________________________________________________________________________
641 Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const {
642 // Checks whether point is inside the EMCal volume, used in AliEMCALv*.cxx
644 // Code uses cylindrical approximation made of inner radius (for speed)
646 // Points behind EMCAl, i.e. R > outer radius, but eta, phi in acceptance
647 // are considered to inside
649 Double_t r=sqrt(x*x+y*y);
651 if ( r > fEnvelop[0] ) {
653 theta = TMath::ATan2(r,z);
658 eta = -TMath::Log(TMath::Tan(theta/2.));
659 if (eta < fArm1EtaMin || eta > fArm1EtaMax)
662 Double_t phi = TMath::ATan2(y,x) * 180./TMath::Pi();
663 if (phi < 0) phi += 360; // phi should go from 0 to 360 in this case
664 if (phi > fArm1PhiMin && phi < fArm1PhiMax)
671 // == Shish-kebab cases ==
673 //________________________________________________________________________________________________
674 Int_t AliEMCALGeometry::GetAbsCellId(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta) const
678 // 13-oct-05; 110 degree case
679 // May 31, 2006; ALICE numbering scheme:
680 // 0 <= nSupMod < fNumberOfSuperModules
681 // 0 <= nModule < fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1)
682 // 0 <= nIphi < fNPHIdiv
683 // 0 <= nIeta < fNETAdiv
684 // 0 <= absid < fNCells
685 static Int_t id=0; // have to change from 0 to fNCells-1
686 if(fKey110DEG == 1 && nSupMod >= 10) { // 110 degree case; last two supermodules
687 id = fNCellsInSupMod*10 + (fNCellsInSupMod/2)*(nSupMod-10);
689 id = fNCellsInSupMod*nSupMod;
691 id += fNCellsInModule *nModule;
692 id += fNPHIdiv *nIphi;
694 if(id<0 || id >= fNCells) {
695 // printf(" wrong numerations !!\n");
696 // printf(" id %6i(will be force to -1)\n", id);
697 // printf(" fNCells %6i\n", fNCells);
698 // printf(" nSupMod %6i\n", nSupMod);
699 // printf(" nModule %6i\n", nModule);
700 // printf(" nIphi %6i\n", nIphi);
701 // printf(" nIeta %6i\n", nIeta);
702 id = -TMath::Abs(id); // if negative something wrong
707 //________________________________________________________________________________________________
708 Bool_t AliEMCALGeometry::CheckAbsCellId(Int_t absId) const
710 // May 31, 2006; only trd1 now
711 if(absId<0 || absId >= fNCells) return kFALSE;
715 //________________________________________________________________________________________________
716 Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nModule,Int_t &nIphi,Int_t &nIeta) const
718 // 21-sep-04; 19-oct-05;
719 // May 31, 2006; ALICE numbering scheme:
722 // absId - cell is as in Geant, 0<= absId < fNCells;
724 // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
725 // nModule - module number in SM, 0<= nModule < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
726 // nIphi - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv;
727 // nIeta - cell number in eta driection inside module; 0<= nIeta < fNETAdiv;
729 static Int_t tmp=0, sm10=0;
730 if(!CheckAbsCellId(absId)) return kFALSE;
732 sm10 = fNCellsInSupMod*10;
733 if(fKey110DEG == 1 && absId >= sm10) { // 110 degree case; last two supermodules
734 nSupMod = (absId-sm10) / (fNCellsInSupMod/2) + 10;
735 tmp = (absId-sm10) % (fNCellsInSupMod/2);
737 nSupMod = absId / fNCellsInSupMod;
738 tmp = absId % fNCellsInSupMod;
741 nModule = tmp / fNCellsInModule;
742 tmp = tmp % fNCellsInModule;
743 nIphi = tmp / fNPHIdiv;
744 nIeta = tmp % fNPHIdiv;
749 //________________________________________________________________________________________________
750 void AliEMCALGeometry::GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, int &iphim, int &ietam) const
752 // added nSupMod; - 19-oct-05 !
753 // Alice numbering scheme - Jun 01,2006
754 // ietam, iphi - indexes of module in two dimensional grid of SM
755 // ietam - have to change from 0 to fNZ-1
756 // iphim - have to change from 0 to nphi-1 (fNPhi-1 or fNPhi/2-1)
759 if(fKey110DEG == 1 && nSupMod>=10) nphi = fNPhi/2;
762 ietam = nModule/nphi;
763 iphim = nModule%nphi;
766 //________________________________________________________________________________________________
767 void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta,
768 int &iphi, int &ieta) const
771 // Added nSupMod; Nov 25, 05
772 // Alice numbering scheme - Jun 01,2006
774 // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
775 // nModule - module number in SM, 0<= nModule < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
776 // nIphi - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv;
777 // nIeta - cell number in eta driection inside module; 0<= nIeta < fNETAdiv;
780 // ieta, iphi - indexes of cell(tower) in two dimensional grid of SM
781 // ieta - have to change from 0 to (fNZ*fNETAdiv-1)
782 // iphi - have to change from 0 to (fNPhi*fNPHIdiv-1 or fNPhi*fNPHIdiv/2-1)
784 static Int_t iphim, ietam;
786 GetModulePhiEtaIndexInSModule(nSupMod,nModule, iphim, ietam);
787 // ieta = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM)
788 ieta = ietam*fNETAdiv + (fNETAdiv - 1 - nIeta); // x(module) = -z(SM)
789 iphi = iphim*fNPHIdiv + nIphi; // y(module) = y(SM)
792 AliDebug(1,Form(" nSupMod %i nModule %i nIphi %i nIeta %i => ieta %i iphi %i\n",
793 nSupMod, nModule, nIphi, nIeta, ieta, iphi));
796 //________________________________________________________________________________________________
797 Int_t AliEMCALGeometry::GetSuperModuleNumber(Int_t absId) const
799 // Return the number of the supermodule given the absolute
800 // ALICE numbering id
802 static Int_t nSupMod, nModule, nIphi, nIeta;
803 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
807 //________________________________________________________________________________________________
808 void AliEMCALGeometry::GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta,
809 Int_t &iphim, Int_t &ietam, Int_t &nModule) const
811 // Transition from cell indexes (ieta,iphi) to module indexes (ietam,iphim, nModule)
813 nphi = GetNumberOfModuleInPhiDirection(nSupMod);
815 ietam = ieta/fNETAdiv;
816 iphim = iphi/fNPHIdiv;
817 nModule = ietam * nphi + iphim;
820 //________________________________________________________________________________________________
821 Int_t AliEMCALGeometry::GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const
823 // Transition from super module number(nSupMod) and cell indexes (ieta,iphi) to absId
824 static Int_t ietam, iphim, nModule;
825 static Int_t nIeta, nIphi; // cell indexes in module
827 GetModuleIndexesFromCellIndexesInSModule(nSupMod, iphi, ieta, ietam, iphim, nModule);
829 nIeta = ieta%fNETAdiv;
830 nIeta = fNETAdiv - 1 - nIeta;
831 nIphi = iphi%fNPHIdiv;
833 return GetAbsCellId(nSupMod, nModule, nIphi, nIeta);
837 // Methods for AliEMCALRecPoint - Feb 19, 2006
838 //________________________________________________________________________________________________
839 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const
841 // Look to see what the relative
842 // position inside a given cell is
844 // Alice numbering scheme - Jun 08, 2006
846 // absId - cell is as in Geant, 0<= absId < fNCells;
848 // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
850 // Shift index taking into account the difference between standard SM
851 // and SM of half size in phi direction
852 const Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
853 static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
854 if(!CheckAbsCellId(absId)) return kFALSE;
856 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
857 GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
859 xr = fCentersOfCellsXDir.At(ieta);
860 zr = fCentersOfCellsEtaDir.At(ieta);
863 yr = fCentersOfCellsPhiDir.At(iphi);
865 yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
867 AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
872 //________________________________________________________________________________________________
873 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t loc[3]) const
875 // Alice numbering scheme - Jun 03, 2006
876 loc[0] = loc[1] = loc[2]=0.0;
877 if(RelPosCellInSModule(absId, loc[0],loc[1],loc[2])) {
883 //________________________________________________________________________________________________
884 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, TVector3 &vloc) const
886 static Double_t loc[3];
887 if(RelPosCellInSModule(absId,loc)) {
888 vloc.SetXYZ(loc[0], loc[1], loc[2]);
894 // Alice numbering scheme - Jun 03, 2006
897 //________________________________________________________________________________________________
898 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t distEff, Double_t &xr, Double_t &yr, Double_t &zr) const
900 // Jul 30, 2007 - taking into account position of shower max
901 // Look to see what the relative
902 // position inside a given cell is
905 // absId - cell is as in Geant, 0<= absId < fNCells;
906 // e - cluster energy
908 // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
910 // Shift index taking into account the difference between standard SM
911 // and SM of half size in phi direction
912 const Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
913 static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
914 static Int_t iphim, ietam;
915 static AliEMCALShishKebabTrd1Module *mod = 0;
917 if(!CheckAbsCellId(absId)) return kFALSE;
919 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
920 GetModulePhiEtaIndexInSModule(nSupMod, nModule, iphim, ietam);
921 GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
923 mod = GetShishKebabModule(ietam);
924 mod->GetPositionAtCenterCellLine(nIeta, distEff, v);
925 xr = v.Y() - fParSM[0];
926 zr = v.X() - fParSM[2];
929 yr = fCentersOfCellsPhiDir.At(iphi);
931 yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
933 AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
938 //________________________________________________________________________________________________
939 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Int_t maxAbsId, Double_t distEff, Double_t &xr, Double_t &yr, Double_t &zr) const
941 // Jul 31, 2007 - taking into account position of shower max and apply coor2.
942 // Look to see what the relative
943 // position inside a given cell is
946 // absId - cell is as in Geant, 0<= absId < fNCells;
947 // maxAbsId - abs id of cell with highest energy
948 // e - cluster energy
950 // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
952 // Shift index taking into account the difference between standard SM
953 // and SM of half size in phi direction
954 const Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
955 static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
956 static Int_t iphim, ietam;
957 static AliEMCALShishKebabTrd1Module *mod = 0;
960 static Int_t nSupModM, nModuleM, nIphiM, nIetaM, iphiM, ietaM;
961 static Int_t iphimM, ietamM, maxAbsIdCopy=-1;
962 static AliEMCALShishKebabTrd1Module *modM = 0;
963 static Double_t distCorr;
965 if(!CheckAbsCellId(absId)) return kFALSE;
967 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
968 GetModulePhiEtaIndexInSModule(nSupMod, nModule, iphim, ietam);
969 GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
970 mod = GetShishKebabModule(ietam);
972 if(absId != maxAbsId) {
974 if(maxAbsIdCopy != maxAbsId) {
975 GetCellIndex(maxAbsId, nSupModM, nModuleM, nIphiM, nIetaM);
976 GetModulePhiEtaIndexInSModule(nSupModM, nModuleM, iphimM, ietamM);
977 GetCellPhiEtaIndexInSModule(nSupModM,nModuleM,nIphiM,nIetaM, iphiM, ietaM);
978 modM = GetShishKebabModule(ietamM); // do I need this ?
979 maxAbsIdCopy = maxAbsId;
982 distCorr = GetEtaModuleSize()*(ietam-ietamM)/TMath::Tan(modM->GetTheta()); // Stay here
983 //printf(" distCorr %f | dist %f | ietam %i -> etamM %i\n", distCorr, dist, ietam, ietamM);
985 // distEff += distCorr;
987 // Bad resolution in this case, strong bias vs phi
989 mod->GetPositionAtCenterCellLine(nIeta, distEff, v); // Stay here
990 xr = v.Y() - fParSM[0];
991 zr = v.X() - fParSM[2];
994 yr = fCentersOfCellsPhiDir.At(iphi);
996 yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
998 AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
1003 //________________________________________________________________________________________________
1004 void AliEMCALGeometry::CreateListOfTrd1Modules()
1006 // Generate the list of Trd1 modules
1007 // which will make up the EMCAL
1010 AliDebug(2,Form(" AliEMCALGeometry::CreateListOfTrd1Modules() started "));
1012 AliEMCALShishKebabTrd1Module *mod=0, *mTmp=0; // current module
1013 if(fShishKebabTrd1Modules == 0) {
1014 fShishKebabTrd1Modules = new TList;
1015 fShishKebabTrd1Modules->SetName("ListOfTRD1");
1016 for(int iz=0; iz< GetNZ(); iz++) {
1018 mod = new AliEMCALShishKebabTrd1Module(TMath::Pi()/2.,this);
1020 mTmp = new AliEMCALShishKebabTrd1Module(*mod);
1023 fShishKebabTrd1Modules->Add(mod);
1026 AliDebug(2,Form(" Already exits : "));
1028 mod = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(fShishKebabTrd1Modules->GetSize()-1);
1029 fEtaMaxOfTRD1 = mod->GetMaxEtaOfModule(0);
1031 AliDebug(2,Form(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n",
1032 fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1));
1034 // Jun 01, 2006 - ALICE numbering scheme
1035 // define grid for cells in eta(z) and x directions in local coordinates system of SM
1036 // Works just for 2x2 case only -- ?? start here
1039 // Define grid for cells in phi(y) direction in local coordinates system of SM
1040 // as for 2X2 as for 3X3 - Nov 8,2006
1042 AliDebug(2,Form(" Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize()));
1043 Int_t ind=0; // this is phi index
1044 Int_t ieta=0, nModule=0, iphiTemp;
1045 Double_t xr=0., zr=0., theta=0., phi=0., eta=0., r=0., x=0.,y=0.;
1047 Double_t ytCenterModule=0.0, ytCenterCell=0.0;
1049 fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
1050 fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
1052 Double_t r0 = GetIPDistance() + GetLongModuleSize()/2.;
1053 for(Int_t it=0; it<fNPhi; it++) { // cycle on modules
1054 ytCenterModule = -fParSM[1] + fPhiModuleSize*(2*it+1)/2; // center of module
1055 for(Int_t ic=0; ic<fNPHIdiv; ic++) { // cycle on cells in module
1057 ytCenterCell = ytCenterModule + fPhiTileSize *(2*ic-1)/2.;
1058 } else if(fNPHIdiv==3){
1059 ytCenterCell = ytCenterModule + fPhiTileSize *(ic-1);
1060 } else if(fNPHIdiv==1){
1061 ytCenterCell = ytCenterModule;
1063 fCentersOfCellsPhiDir.AddAt(ytCenterCell,ind);
1064 // Define grid on phi direction
1065 // Grid is not the same for different eta bin;
1066 // Effect is small but is still here
1067 phi = TMath::ATan2(ytCenterCell, r0);
1068 fPhiCentersOfCells.AddAt(phi, ind);
1070 AliDebug(2,Form(" ind %2.2i : y %8.3f ", ind, fCentersOfCellsPhiDir.At(ind)));
1075 fCentersOfCellsEtaDir.Set(fNZ *fNETAdiv);
1076 fCentersOfCellsXDir.Set(fNZ *fNETAdiv);
1077 fEtaCentersOfCells.Set(fNZ *fNETAdiv * fNPhi*fNPHIdiv);
1078 AliDebug(2,Form(" Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize()));
1079 for(Int_t it=0; it<fNZ; it++) {
1080 AliEMCALShishKebabTrd1Module *trd1 = GetShishKebabModule(it);
1082 for(Int_t ic=0; ic<fNETAdiv; ic++) {
1084 trd1->GetCenterOfCellInLocalCoordinateofSM(ic, xr, zr); // case of 2X2
1085 GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1087 trd1->GetCenterOfCellInLocalCoordinateofSM_3X3(ic, xr, zr); // case of 3X3
1088 GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1090 trd1->GetCenterOfCellInLocalCoordinateofSM_1X1(xr, zr); // case of 1X1
1091 GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1093 fCentersOfCellsXDir.AddAt(float(xr) - fParSM[0],ieta);
1094 fCentersOfCellsEtaDir.AddAt(float(zr) - fParSM[2],ieta);
1095 // Define grid on eta direction for each bin in phi
1096 for(int iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
1097 x = xr + trd1->GetRadius();
1098 y = fCentersOfCellsPhiDir[iphi];
1099 r = TMath::Sqrt(x*x + y*y + zr*zr);
1100 theta = TMath::ACos(zr/r);
1101 eta = AliEMCALShishKebabTrd1Module::ThetaToEta(theta);
1102 // ind = ieta*fCentersOfCellsPhiDir.GetSize() + iphi;
1103 ind = iphi*fCentersOfCellsEtaDir.GetSize() + ieta;
1104 fEtaCentersOfCells.AddAt(eta, ind);
1106 //printf(" ieta %i : xr + trd1->GetRadius() %f : zr %f : eta %f \n", ieta, xr + trd1->GetRadius(), zr, eta);
1109 for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1110 AliDebug(2,Form(" ind %2.2i : z %8.3f : x %8.3f", i+1,
1111 fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i)));
1116 //________________________________________________________________________________________________
1117 void AliEMCALGeometry::GetGlobal(const Double_t *loc, Double_t *glob, int ind) const
1119 // Figure out the global numbering
1120 // of a given supermodule from the
1121 // local numbering and the transformation
1122 // matrix stored by the geometry manager (allows for misaligned
1125 if(ind>=0 && ind < GetNumberOfSuperModules()) {
1126 TString volpath = "ALIC_1/XEN1_1/SMOD_";
1129 if(GetKey110DEG() && ind>=10) {
1130 volpath = "ALIC_1/XEN1_1/SM10_";
1131 volpath += ind-10+1;
1134 if(!gGeoManager->cd(volpath.Data()))
1135 AliFatal(Form("AliEMCALGeometry::GeoManager cannot find path %s!",volpath.Data()));
1137 TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1139 m->LocalToMaster(loc, glob);
1141 AliFatal("Geo matrixes are not loaded \n") ;
1146 //________________________________________________________________________________________________
1147 void AliEMCALGeometry::GetGlobal(const TVector3 &vloc, TVector3 &vglob, int ind) const
1149 //Figure out the global numbering
1150 //of a given supermodule from the
1151 //local numbering given a 3-vector location
1153 static Double_t tglob[3], tloc[3];
1155 GetGlobal(tloc, tglob, ind);
1156 vglob.SetXYZ(tglob[0], tglob[1], tglob[2]);
1159 //________________________________________________________________________________________________
1160 void AliEMCALGeometry::GetGlobal(Int_t absId , double glob[3]) const
1162 // Alice numbering scheme - Jun 03, 2006
1163 static Int_t nSupMod, nModule, nIphi, nIeta;
1164 static double loc[3];
1166 if (!gGeoManager || !gGeoManager->IsClosed()) {
1167 AliError("Can't get the global coordinates! gGeoManager doesn't exist or it is still open!");
1171 glob[0]=glob[1]=glob[2]=0.0; // bad case
1172 if(RelPosCellInSModule(absId, loc)) {
1173 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1175 TString volpath = "ALIC_1/XEN1_1/SMOD_";
1176 volpath += (nSupMod+1);
1178 if(GetKey110DEG() && nSupMod>=10) {
1179 volpath = "ALIC_1/XEN1_1/SM10_";
1180 volpath += (nSupMod-10+1);
1182 if(!gGeoManager->cd(volpath.Data()))
1183 AliFatal(Form("GeoManager cannot find path %s!",volpath.Data()));
1185 TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1187 m->LocalToMaster(loc, glob);
1189 AliFatal("Geo matrixes are not loaded \n") ;
1194 //___________________________________________________________________
1195 void AliEMCALGeometry::GetGlobal(Int_t absId , TVector3 &vglob) const
1197 // Alice numbering scheme - Jun 03, 2006
1198 static Double_t glob[3];
1200 GetGlobal(absId, glob);
1201 vglob.SetXYZ(glob[0], glob[1], glob[2]);
1205 //____________________________________________________________________________
1206 void AliEMCALGeometry::GetGlobal(const AliRecPoint* /*rp*/, TVector3& /* vglob */) const
1208 AliFatal(Form("Please use GetGlobalEMCAL(recPoint,gpos) instead of GetGlobal!"));
1211 //_________________________________________________________________________________
1212 void AliEMCALGeometry::GetGlobalEMCAL(const AliEMCALRecPoint *rp, TVector3 &vglob) const
1214 // Figure out the global numbering
1215 // of a given supermodule from the
1216 // local numbering for RecPoints
1218 static TVector3 vloc;
1219 static Int_t nSupMod, nModule, nIphi, nIeta;
1221 const AliEMCALRecPoint *rpTmp = rp;
1222 const AliEMCALRecPoint *rpEmc = rpTmp;
1224 GetCellIndex(rpEmc->GetAbsId(0), nSupMod, nModule, nIphi, nIeta);
1225 rpTmp->GetLocalPosition(vloc);
1226 GetGlobal(vloc, vglob, nSupMod);
1229 //________________________________________________________________________________________________
1230 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Double_t &eta,Double_t &phi) const
1232 // Nov 16, 2006- float to double
1233 // version for TRD1 only
1234 static TVector3 vglob;
1235 GetGlobal(absId, vglob);
1240 //________________________________________________________________________________________________
1241 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Float_t &eta,Float_t &phi) const
1243 // Nov 16,2006 - should be discard in future
1244 static TVector3 vglob;
1245 GetGlobal(absId, vglob);
1246 eta = float(vglob.Eta());
1247 phi = float(vglob.Phi());
1250 //________________________________________________________________________________________________
1251 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSM(Int_t nSupMod, Double_t &phiMin, Double_t &phiMax) const
1253 // 0<= nSupMod <=11; phi in rad
1255 if(nSupMod<0 || nSupMod >11) return kFALSE;
1257 phiMin = fPhiBoundariesOfSM[2*i];
1258 phiMax = fPhiBoundariesOfSM[2*i+1];
1262 //________________________________________________________________________________________________
1263 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const
1265 // 0<= nPhiSec <=4; phi in rad
1266 // 0; gap boundaries between 0th&2th | 1th&3th SM
1267 // 1; gap boundaries between 2th&4th | 3th&5th SM
1268 // 2; gap boundaries between 4th&6th | 5th&7th SM
1269 // 3; gap boundaries between 6th&8th | 7th&9th SM
1270 // 4; gap boundaries between 8th&10th | 9th&11th SM
1271 if(nPhiSec<0 || nPhiSec >4) return kFALSE;
1272 phiMin = fPhiBoundariesOfSM[2*nPhiSec+1];
1273 phiMax = fPhiBoundariesOfSM[2*nPhiSec+2];
1277 //________________________________________________________________________________________________
1278 Bool_t AliEMCALGeometry::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const
1280 // Return false if phi belongs a phi cracks between SM
1284 if(TMath::Abs(eta) > fEtaMaxOfTRD1) return kFALSE;
1286 phi = TVector2::Phi_0_2pi(phi); // move phi to (0,2pi) boundaries
1287 for(i=0; i<6; i++) {
1288 if(phi>=fPhiBoundariesOfSM[2*i] && phi<=fPhiBoundariesOfSM[2*i+1]) {
1290 if(eta < 0.0) nSupMod++;
1291 AliDebug(1,Form("eta %f phi %f(%5.2f) : nSupMod %i : #bound %i", eta,phi,phi*TMath::RadToDeg(), nSupMod,i));
1298 //________________________________________________________________________________________________
1299 Bool_t AliEMCALGeometry::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi, Int_t &absId) const
1302 // stay here - phi problem as usual
1303 static Int_t nSupMod, i, ieta, iphi, etaShift, nphi;
1304 static Double_t absEta=0.0, d=0.0, dmin=0.0, phiLoc;
1305 absId = nSupMod = - 1;
1306 if(SuperModuleNumberFromEtaPhi(eta, phi, nSupMod)) {
1308 phi = TVector2::Phi_0_2pi(phi);
1309 phiLoc = phi - fPhiCentersOfSM[nSupMod/2];
1310 nphi = fPhiCentersOfCells.GetSize();
1312 phiLoc = phi - 190.*TMath::DegToRad();
1316 dmin = TMath::Abs(fPhiCentersOfCells[0]-phiLoc);
1318 for(i=1; i<nphi; i++) {
1319 d = TMath::Abs(fPhiCentersOfCells[i] - phiLoc);
1324 // printf(" i %i : d %f : dmin %f : fPhiCentersOfCells[i] %f \n", i, d, dmin, fPhiCentersOfCells[i]);
1326 // odd SM are turned with respect of even SM - reverse indexes
1327 AliDebug(2,Form(" iphi %i : dmin %f (phi %f, phiLoc %f ) ", iphi, dmin, phi, phiLoc));
1329 absEta = TMath::Abs(eta);
1330 etaShift = iphi*fCentersOfCellsEtaDir.GetSize();
1331 dmin = TMath::Abs(fEtaCentersOfCells[etaShift]-absEta);
1333 for(i=1; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1334 d = TMath::Abs(fEtaCentersOfCells[i+etaShift] - absEta);
1340 AliDebug(2,Form(" ieta %i : dmin %f (eta=%f) : nSupMod %i ", ieta, dmin, eta, nSupMod));
1342 if(eta<0) iphi = (nphi-1) - iphi;
1343 absId = GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);
1350 //________________________________________________________________________________________________
1351 AliEMCALShishKebabTrd1Module* AliEMCALGeometry::GetShishKebabModule(Int_t neta) const
1353 //This method was too long to be
1354 //included in the header file - the
1355 //rule checker complained about it's
1356 //length, so we move it here. It returns the
1357 //shishkebabmodule at a given eta index point.
1359 static AliEMCALShishKebabTrd1Module* trd1=0;
1360 if(fShishKebabTrd1Modules && neta>=0 && neta<fShishKebabTrd1Modules->GetSize()) {
1361 trd1 = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(neta);
1366 //________________________________________________________________________________________________
1367 Int_t AliEMCALGeometry::GetAbsTRUNumberFromNumberInSm(const Int_t row, const Int_t col, const Int_t sm)
1369 Int_t itru = row + col*GetNModulesInTRUPhi() + sm*GetNTRU();
1370 // printf(" GetAbsTRUNumberFromNumberInSm : row %2i col %2i sm %2i -> itru %2i\n", row, col, sm, itru);
1374 //________________________________________________________________________________________________
1375 void AliEMCALGeometry::Browse(TBrowser* b)
1377 //Browse the modules
1378 if(fShishKebabTrd1Modules) b->Add(fShishKebabTrd1Modules);
1381 //________________________________________________________________________________________________
1382 Bool_t AliEMCALGeometry::IsFolder() const
1384 //Check if fShishKebabTrd1Modules is in folder
1385 if(fShishKebabTrd1Modules) return kTRUE;
1389 //________________________________________________________________________________________________
1390 Double_t AliEMCALGeometry::GetPhiCenterOfSM(Int_t nsupmod) const
1392 //returns center of supermodule in phi
1394 return fPhiCentersOfSM[i];
1397 //____________________________________________________________________________
1398 Bool_t AliEMCALGeometry::Impact(const TParticle * particle) const
1400 // Tells if a particle enters EMCAL
1403 TVector3 vtx(particle->Vx(),particle->Vy(),particle->Vz());
1404 TVector3 vimpact(0,0,0);
1405 ImpactOnEmcal(vtx,particle->Theta(),particle->Phi(),AbsID,vimpact);
1410 //____________________________________________________________________________
1411 void AliEMCALGeometry::ImpactOnEmcal(TVector3 vtx, Double_t theta, Double_t phi,
1412 Int_t & absId, TVector3 & vimpact) const
1414 // calculates the impact coordinates on EMCAL (centre of a tower/not on EMCAL surface)
1415 // of a neutral particle
1416 // emitted in the vertex vtx[3] with direction theta and phi in the ALICE global coordinate system
1418 TVector3 p(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) ;
1420 vimpact.SetXYZ(0,0,0);
1422 if(phi==0 || theta==0) return;
1425 Double_t factor = (GetIPDistance()-vtx[1])/p[1];
1426 direction = vtx + factor*p;
1429 AliFatal("Geo manager not initialized\n");
1431 //from particle direction -> tower hitted
1432 GetAbsCellIdFromEtaPhi(direction.Eta(),direction.Phi(),absId);
1434 //tower absID hitted -> tower/module plane (evaluated at the center of the tower)
1435 Int_t nSupMod, nModule, nIphi, nIeta;
1436 Double_t loc[3],loc2[3],loc3[3];
1437 Double_t glob[3]={},glob2[3]={},glob3[3]={};
1439 if(!RelPosCellInSModule(absId,loc)) return;
1441 //loc is cell center of tower
1442 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1444 //look at 2 neighbours-s cell using nIphi={0,1} and nIeta={0,1}
1445 Int_t nIphi2,nIeta2,absId2,absId3;
1446 if(nIeta==0) nIeta2=1;
1448 absId2=GetAbsCellId(nSupMod,nModule,nIphi,nIeta2);
1449 if(nIphi==0) nIphi2=1;
1451 absId3=GetAbsCellId(nSupMod,nModule,nIphi2,nIeta);
1453 //2nd point on emcal cell plane
1454 if(!RelPosCellInSModule(absId2,loc2)) return;
1456 //3rd point on emcal cell plane
1457 if(!RelPosCellInSModule(absId3,loc3)) return;
1459 TString volpath = "ALIC_1/XEN1_1/SMOD_";
1460 volpath += (nSupMod+1);
1462 if(GetKey110DEG() && nSupMod>=10) {
1463 volpath = "ALIC_1/XEN1_1/SM10_";
1464 volpath += (nSupMod-10+1);
1466 if(!gGeoManager->cd(volpath.Data())){
1467 AliFatal(Form("GeoManager cannot find path %s!",volpath.Data()))
1470 TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1472 m->LocalToMaster(loc, glob);
1473 m->LocalToMaster(loc2, glob2);
1474 m->LocalToMaster(loc3, glob3);
1476 AliFatal("Geo matrixes are not loaded \n") ;
1479 //Equation of Plane from glob,glob2,glob3 (Ax+By+Cz+D=0)
1480 Double_t A = glob[1]*(glob2[2]-glob3[2]) + glob2[1]*(glob3[2]-glob[2]) + glob3[1]*(glob[2]-glob2[2]);
1481 Double_t B = glob[2]*(glob2[0]-glob3[0]) + glob2[2]*(glob3[0]-glob[0]) + glob3[2]*(glob[0]-glob2[0]);
1482 Double_t C = glob[0]*(glob2[1]-glob3[1]) + glob2[0]*(glob3[1]-glob[1]) + glob3[0]*(glob[1]-glob2[1]);
1483 Double_t D = glob[0]*(glob2[1]*glob3[2]-glob3[1]*glob2[2]) + glob2[0]*(glob3[1]*glob[2]-glob[1]*glob3[2]) + glob3[0]*(glob[1]*glob2[2]-glob2[1]*glob[2]);
1486 //shift equation of plane from tower/module center to surface along vector (A,B,C) normal to tower/module plane
1487 Double_t dist = GetLongModuleSize()/2.;
1488 Double_t norm = TMath::Sqrt(A*A+B*B+C*C);
1489 Double_t glob4[3]={};
1490 TVector3 dir(A,B,C);
1491 TVector3 point(glob[0],glob[1],glob[2]);
1492 if(point.Dot(dir)<0) dist*=-1;
1493 glob4[0]=glob[0]-dist*A/norm;
1494 glob4[1]=glob[1]-dist*B/norm;
1495 glob4[2]=glob[2]-dist*C/norm;
1496 D = glob4[0]*A + glob4[1]*B + glob4[2]*C ;
1499 //Line determination (2 points for equation of line : vtx and direction)
1500 //impact between line (particle) and plane (module/tower plane)
1501 Double_t den = A*(vtx(0)-direction(0)) + B*(vtx(1)-direction(1)) + C*(vtx(2)-direction(2));
1503 printf("ImpactOnEmcal() No solution :\n");
1507 Double_t length = A*vtx(0)+B*vtx(1)+C*vtx(2)+D;
1510 vimpact.SetXYZ(vtx(0)+length*(direction(0)-vtx(0)),vtx(1)+length*(direction(1)-vtx(1)),vtx(2)+length*(direction(2)-vtx(2)));
1512 //shift vimpact from tower/module surface to center along vector (A,B,C) normal to tower/module plane
1513 vimpact.SetXYZ(vimpact(0)+dist*A/norm,vimpact(1)+dist*B/norm,vimpact(2)+dist*C/norm);