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 void AliEMCALGeometry::PrintGeometry()
401 // Separate routine is callable from broswer; Nov 7,2006
402 printf("\nInit: geometry of EMCAL named %s :\n", fGeoName.Data());
404 for(Int_t i=0; i<fArrayOpts->GetEntries(); i++){
405 TObjString *o = (TObjString*)fArrayOpts->At(i);
406 printf(" %i : %s \n", i, o->String().Data());
409 printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
410 printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f -> for EMCAL envelope only\n",
411 GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() );
413 printf( " ECAL : %d x (%f cm Pb, %f cm Sc) \n",
414 GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ;
415 printf(" fSampling %5.2f \n", fSampling );
416 printf(" fIPDistance %6.3f cm \n", fIPDistance);
417 printf(" fNPhi %i | fNZ %i \n", fNPhi, fNZ);
418 printf(" fNCellsInModule %i : fNCellsInSupMod %i : fNCells %i\n",fNCellsInModule, fNCellsInSupMod, fNCells);
419 printf(" X:Y module size %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize);
420 printf(" X:Y tile size %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize);
421 printf(" #of sampling layers %i(fNECLayers) \n", fNECLayers);
422 printf(" fLongModuleSize %6.3f cm \n", fLongModuleSize);
423 printf(" #supermodule in phi direction %i \n", fNPhiSuperModule );
424 printf(" fILOSS %i : fIHADR %i \n", fILOSS, fIHADR);
425 printf(" fTrd1Angle %7.4f\n", fTrd1Angle);
426 printf(" f2Trd1Dx2 %7.4f\n", f2Trd1Dx2);
427 printf("SM dimensions(TRD1) : dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n",
428 fParSM[0],fParSM[1],fParSM[2]);
429 printf(" fPhiGapForSM %7.4f cm (%7.4f <- phi size in degree)\n",
430 fPhiGapForSM, TMath::ATan2(fPhiGapForSM,fIPDistance)*TMath::RadToDeg());
431 if(GetKey110DEG()) printf(" Last two modules have size 10 degree in phi (180<phi<190)\n");
432 printf(" phi SM boundaries \n");
433 for(int i=0; i<fPhiBoundariesOfSM.GetSize()/2.; i++) {
434 printf(" %i : %7.5f(%7.2f) -> %7.5f(%7.2f) : center %7.5f(%7.2f) \n", i,
435 fPhiBoundariesOfSM[2*i], fPhiBoundariesOfSM[2*i]*TMath::RadToDeg(),
436 fPhiBoundariesOfSM[2*i+1], fPhiBoundariesOfSM[2*i+1]*TMath::RadToDeg(),
437 fPhiCentersOfSM[i], fPhiCentersOfSM[i]*TMath::RadToDeg());
439 printf(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n",
440 fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1);
442 printf("\n Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize());
443 for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
444 printf(" ind %2.2i : z %8.3f : x %8.3f \n", i,
445 fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i));
446 int ind=0; // Nov 21,2006
447 for(Int_t iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
448 ind = iphi*fCentersOfCellsEtaDir.GetSize() + i;
449 printf("%6.4f ", fEtaCentersOfCells[ind]);
450 if((iphi+1)%12 == 0) printf("\n");
456 printf("\n Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize());
457 for(Int_t i=0; i<fCentersOfCellsPhiDir.GetSize(); i++) {
458 double phi=fPhiCentersOfCells.At(i);
459 printf(" ind %2.2i : y %8.3f : phi %7.5f(%6.2f) \n", i, fCentersOfCellsPhiDir.At(i),
460 phi, phi*TMath::RadToDeg());
465 //______________________________________________________________________
466 void AliEMCALGeometry::PrintCellIndexes(Int_t absId, int pri, char *tit)
469 Int_t nSupMod, nModule, nIphi, nIeta;
473 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
474 printf(" %s | absId : %i -> nSupMod %i nModule %i nIphi %i nIeta %i \n", tit, absId, nSupMod, nModule, nIphi, nIeta);
476 GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
477 printf(" local SM index : iphi %i : ieta %i \n", iphi,ieta);
478 GetGlobal(absId, vg);
479 printf(" vglob : mag %7.2f : perp %7.2f : z %7.2f : eta %6.4f : phi %6.4f(%6.2f) \n",
480 vg.Mag(), vg.Perp(), vg.Z(), vg.Eta(), vg.Phi(), vg.Phi()*TMath::RadToDeg());
484 //______________________________________________________________________
485 void AliEMCALGeometry::CheckAdditionalOptions()
488 // Additional options that
489 // can be used to select
490 // the specific geometry of
493 // adeed allILOSS= and allIHADR= for MIP investigation
494 fArrayOpts = new TObjArray;
495 Int_t nopt = AliEMCALHistoUtilities::ParseString(fGeoName, *fArrayOpts);
496 if(nopt==1) { // no aditional option(s)
497 fArrayOpts->Delete();
502 for(Int_t i=1; i<nopt; i++){
503 TObjString *o = (TObjString*)fArrayOpts->At(i);
505 TString addOpt = o->String();
507 for(Int_t j=0; j<fNAdditionalOpts; j++) {
508 TString opt = fAdditionalOpts[j];
509 if(addOpt.Contains(opt,TString::kIgnoreCase)) {
515 AliDebug(2,Form("<E> option |%s| unavailable : ** look to the file AliEMCALGeometry.h **\n",
519 AliDebug(2,Form("<I> option |%s| is valid : number %i : |%s|\n",
520 addOpt.Data(), indj, fAdditionalOpts[indj]));
521 if (addOpt.Contains("NL=",TString::kIgnoreCase)) {// number of sampling layers
522 sscanf(addOpt.Data(),"NL=%i", &fNECLayers);
523 AliDebug(2,Form(" fNECLayers %i (new) \n", fNECLayers));
524 } else if(addOpt.Contains("PBTH=",TString::kIgnoreCase)) {//Thickness of the Pb(fECPbRadThicknes)
525 sscanf(addOpt.Data(),"PBTH=%f", &fECPbRadThickness);
526 } else if(addOpt.Contains("SCTH=",TString::kIgnoreCase)) {//Thickness of the Sc(fECScintThick)
527 sscanf(addOpt.Data(),"SCTH=%f", &fECScintThick);
528 } else if(addOpt.Contains("LATSS=",TString::kIgnoreCase)) {// Thickness of lateral steel strip (fLateralSteelStrip)
529 sscanf(addOpt.Data(),"LATSS=%f", &fLateralSteelStrip);
530 AliDebug(2,Form(" fLateralSteelStrip %f (new) \n", fLateralSteelStrip));
531 } else if(addOpt.Contains("ILOSS=",TString::kIgnoreCase)) {// As in Geant
532 sscanf(addOpt.Data(),"ALLILOSS=%i", &fILOSS);
533 AliDebug(2,Form(" fILOSS %i \n", fILOSS));
534 } else if(addOpt.Contains("IHADR=",TString::kIgnoreCase)) {// As in Geant
535 sscanf(addOpt.Data(),"ALLIHADR=%i", &fIHADR);
536 AliDebug(2,Form(" fIHADR %i \n", fIHADR));
542 void AliEMCALGeometry::DefineSamplingFraction()
545 // Look http://rhic.physics.wayne.edu/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
546 // Keep for compatibilty
548 if(fNECLayers == 69) { // 10% layer reduction
550 } else if(fNECLayers == 61) { // 20% layer reduction
552 } else if(fNECLayers == 77) {
553 if (fECScintThick>0.159 && fECScintThick<0.161) { // original sampling fraction, equal layers
554 fSampling = 12.327; // fECScintThick = fECPbRadThickness = 0.160;
555 } else if (fECScintThick>0.175 && fECScintThick<0.177) { // 10% Pb thicknes reduction
556 fSampling = 10.5; // fECScintThick = 0.176, fECPbRadThickness=0.144;
557 } else if(fECScintThick>0.191 && fECScintThick<0.193) { // 20% Pb thicknes reduction
558 fSampling = 8.93; // fECScintThick = 0.192, fECPbRadThickness=0.128;
564 //______________________________________________________________________
565 void AliEMCALGeometry::GetModulePhiEtaIndexInSModuleFromTRUIndex(const Int_t itru, const Int_t iphitru, const Int_t ietatru, Int_t &iphiSM, Int_t &ietaSM) const
568 // This method transforms the (eta,phi) index of module in a
569 // TRU matrix into Super Module (eta,phi) index.
571 // Calculate in which row and column where the TRU are
574 Int_t col = itru/ fNTRUPhi ; // indexes of TRU in SM
575 Int_t row = itru - col*fNTRUPhi ;
577 iphiSM = fNModulesInTRUPhi*row + iphitru ;
578 ietaSM = fNModulesInTRUEta*col + ietatru ;
579 //printf(" GetModulePhiEtaIndexInSModuleFromTRUIndex : itru %2i iphitru %2i ietatru %2i iphiSM %2i ietaSM %2i \n",
580 // itru, iphitru, ietatru, iphiSM, ietaSM);
583 //______________________________________________________________________
584 AliEMCALGeometry * AliEMCALGeometry::GetInstance(){
585 // Returns the pointer of the unique instance
587 AliEMCALGeometry * rv = static_cast<AliEMCALGeometry *>( fgGeom );
591 //______________________________________________________________________
592 AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
593 const Text_t* title){
594 // Returns the pointer of the unique instance
596 AliEMCALGeometry * rv = 0;
598 if ( strcmp(name,"") == 0 ) { // get default geometry
599 fgGeom = new AliEMCALGeometry(fgDefaultGeometryName, title);
601 fgGeom = new AliEMCALGeometry(name, title);
602 } // end if strcmp(name,"")
603 if ( fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
610 if ( strcmp(fgGeom->GetName(), name) != 0) {
611 printf("\ncurrent geometry is %s : ", fgGeom->GetName());
612 printf(" you cannot call %s ",name);
614 rv = (AliEMCALGeometry *) fgGeom;
620 //_____________________________________________________________________________
621 Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const {
622 // Checks whether point is inside the EMCal volume, used in AliEMCALv*.cxx
624 // Code uses cylindrical approximation made of inner radius (for speed)
626 // Points behind EMCAl, i.e. R > outer radius, but eta, phi in acceptance
627 // are considered to inside
629 Double_t r=sqrt(x*x+y*y);
631 if ( r > fEnvelop[0] ) {
633 theta = TMath::ATan2(r,z);
638 eta = -TMath::Log(TMath::Tan(theta/2.));
639 if (eta < fArm1EtaMin || eta > fArm1EtaMax)
642 Double_t phi = TMath::ATan2(y,x) * 180./TMath::Pi();
643 if (phi < 0) phi += 360; // phi should go from 0 to 360 in this case
644 if (phi > fArm1PhiMin && phi < fArm1PhiMax)
651 // == Shish-kebab cases ==
653 //________________________________________________________________________________________________
654 Int_t AliEMCALGeometry::GetAbsCellId(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta) const
658 // 13-oct-05; 110 degree case
659 // May 31, 2006; ALICE numbering scheme:
660 // 0 <= nSupMod < fNumberOfSuperModules
661 // 0 <= nModule < fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1)
662 // 0 <= nIphi < fNPHIdiv
663 // 0 <= nIeta < fNETAdiv
664 // 0 <= absid < fNCells
665 static Int_t id=0; // have to change from 0 to fNCells-1
666 if(fKey110DEG == 1 && nSupMod >= 10) { // 110 degree case; last two supermodules
667 id = fNCellsInSupMod*10 + (fNCellsInSupMod/2)*(nSupMod-10);
669 id = fNCellsInSupMod*nSupMod;
671 id += fNCellsInModule *nModule;
672 id += fNPHIdiv *nIphi;
674 if(id<0 || id >= fNCells) {
675 // printf(" wrong numerations !!\n");
676 // printf(" id %6i(will be force to -1)\n", id);
677 // printf(" fNCells %6i\n", fNCells);
678 // printf(" nSupMod %6i\n", nSupMod);
679 // printf(" nModule %6i\n", nModule);
680 // printf(" nIphi %6i\n", nIphi);
681 // printf(" nIeta %6i\n", nIeta);
682 id = -TMath::Abs(id); // if negative something wrong
687 //________________________________________________________________________________________________
688 Bool_t AliEMCALGeometry::CheckAbsCellId(Int_t absId) const
690 // May 31, 2006; only trd1 now
691 if(absId<0 || absId >= fNCells) return kFALSE;
695 //________________________________________________________________________________________________
696 Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nModule,Int_t &nIphi,Int_t &nIeta) const
698 // 21-sep-04; 19-oct-05;
699 // May 31, 2006; ALICE numbering scheme:
702 // absId - cell is as in Geant, 0<= absId < fNCells;
704 // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
705 // nModule - module number in SM, 0<= nModule < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
706 // nIphi - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv;
707 // nIeta - cell number in eta driection inside module; 0<= nIeta < fNETAdiv;
709 static Int_t tmp=0, sm10=0;
710 if(!CheckAbsCellId(absId)) return kFALSE;
712 sm10 = fNCellsInSupMod*10;
713 if(fKey110DEG == 1 && absId >= sm10) { // 110 degree case; last two supermodules
714 nSupMod = (absId-sm10) / (fNCellsInSupMod/2) + 10;
715 tmp = (absId-sm10) % (fNCellsInSupMod/2);
717 nSupMod = absId / fNCellsInSupMod;
718 tmp = absId % fNCellsInSupMod;
721 nModule = tmp / fNCellsInModule;
722 tmp = tmp % fNCellsInModule;
723 nIphi = tmp / fNPHIdiv;
724 nIeta = tmp % fNPHIdiv;
729 //________________________________________________________________________________________________
730 void AliEMCALGeometry::GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, int &iphim, int &ietam) const
732 // added nSupMod; - 19-oct-05 !
733 // Alice numbering scheme - Jun 01,2006
734 // ietam, iphi - indexes of module in two dimensional grid of SM
735 // ietam - have to change from 0 to fNZ-1
736 // iphim - have to change from 0 to nphi-1 (fNPhi-1 or fNPhi/2-1)
739 if(fKey110DEG == 1 && nSupMod>=10) nphi = fNPhi/2;
742 ietam = nModule/nphi;
743 iphim = nModule%nphi;
746 //________________________________________________________________________________________________
747 void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta,
748 int &iphi, int &ieta) const
751 // Added nSupMod; Nov 25, 05
752 // Alice numbering scheme - Jun 01,2006
754 // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
755 // nModule - module number in SM, 0<= nModule < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
756 // nIphi - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv;
757 // nIeta - cell number in eta driection inside module; 0<= nIeta < fNETAdiv;
760 // ieta, iphi - indexes of cell(tower) in two dimensional grid of SM
761 // ieta - have to change from 0 to (fNZ*fNETAdiv-1)
762 // iphi - have to change from 0 to (fNPhi*fNPHIdiv-1 or fNPhi*fNPHIdiv/2-1)
764 static Int_t iphim, ietam;
766 GetModulePhiEtaIndexInSModule(nSupMod,nModule, iphim, ietam);
767 // ieta = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM)
768 ieta = ietam*fNETAdiv + (fNETAdiv - 1 - nIeta); // x(module) = -z(SM)
769 iphi = iphim*fNPHIdiv + nIphi; // y(module) = y(SM)
772 AliDebug(1,Form(" nSupMod %i nModule %i nIphi %i nIeta %i => ieta %i iphi %i\n",
773 nSupMod, nModule, nIphi, nIeta, ieta, iphi));
776 //________________________________________________________________________________________________
777 Int_t AliEMCALGeometry::GetSuperModuleNumber(Int_t absId) const
779 // Return the number of the supermodule given the absolute
780 // ALICE numbering id
782 static Int_t nSupMod, nModule, nIphi, nIeta;
783 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
787 //________________________________________________________________________________________________
788 void AliEMCALGeometry::GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta,
789 Int_t &iphim, Int_t &ietam, Int_t &nModule) const
791 // Transition from cell indexes (ieta,iphi) to module indexes (ietam,iphim, nModule)
793 nphi = GetNumberOfModuleInPhiDirection(nSupMod);
795 ietam = ieta/fNETAdiv;
796 iphim = iphi/fNPHIdiv;
797 nModule = ietam * nphi + iphim;
800 //________________________________________________________________________________________________
801 Int_t AliEMCALGeometry::GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const
803 // Transition from super module number(nSupMod) and cell indexes (ieta,iphi) to absId
804 static Int_t ietam, iphim, nModule;
805 static Int_t nIeta, nIphi; // cell indexes in module
807 GetModuleIndexesFromCellIndexesInSModule(nSupMod, iphi, ieta, ietam, iphim, nModule);
809 nIeta = ieta%fNETAdiv;
810 nIeta = fNETAdiv - 1 - nIeta;
811 nIphi = iphi%fNPHIdiv;
813 return GetAbsCellId(nSupMod, nModule, nIphi, nIeta);
817 // Methods for AliEMCALRecPoint - Feb 19, 2006
818 //________________________________________________________________________________________________
819 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const
821 // Look to see what the relative
822 // position inside a given cell is
824 // Alice numbering scheme - Jun 08, 2006
826 // absId - cell is as in Geant, 0<= absId < fNCells;
828 // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
830 // Shift index taking into account the difference between standard SM
831 // and SM of half size in phi direction
832 const Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
833 static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
834 if(!CheckAbsCellId(absId)) return kFALSE;
836 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
837 GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
839 xr = fCentersOfCellsXDir.At(ieta);
840 zr = fCentersOfCellsEtaDir.At(ieta);
843 yr = fCentersOfCellsPhiDir.At(iphi);
845 yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
847 AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
852 //________________________________________________________________________________________________
853 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t loc[3]) const
855 // Alice numbering scheme - Jun 03, 2006
856 loc[0] = loc[1] = loc[2]=0.0;
857 if(RelPosCellInSModule(absId, loc[0],loc[1],loc[2])) {
863 //________________________________________________________________________________________________
864 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, TVector3 &vloc) const
866 static Double_t loc[3];
867 if(RelPosCellInSModule(absId,loc)) {
868 vloc.SetXYZ(loc[0], loc[1], loc[2]);
874 // Alice numbering scheme - Jun 03, 2006
877 //________________________________________________________________________________________________
878 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t distEff, Double_t &xr, Double_t &yr, Double_t &zr) const
880 // Jul 30, 2007 - taking into account position of shower max
881 // Look to see what the relative
882 // position inside a given cell is
885 // absId - cell is as in Geant, 0<= absId < fNCells;
886 // e - cluster energy
888 // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
890 // Shift index taking into account the difference between standard SM
891 // and SM of half size in phi direction
892 const Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
893 static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
894 static Int_t iphim, ietam;
895 static AliEMCALShishKebabTrd1Module *mod = 0;
897 if(!CheckAbsCellId(absId)) return kFALSE;
899 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
900 GetModulePhiEtaIndexInSModule(nSupMod, nModule, iphim, ietam);
901 GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
903 mod = GetShishKebabModule(ietam);
904 mod->GetPositionAtCenterCellLine(nIeta, distEff, v);
905 xr = v.Y() - fParSM[0];
906 zr = v.X() - fParSM[2];
909 yr = fCentersOfCellsPhiDir.At(iphi);
911 yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
913 AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
918 //________________________________________________________________________________________________
919 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Int_t maxAbsId, Double_t distEff, Double_t &xr, Double_t &yr, Double_t &zr) const
921 // Jul 31, 2007 - taking into account position of shower max and apply coor2.
922 // Look to see what the relative
923 // position inside a given cell is
926 // absId - cell is as in Geant, 0<= absId < fNCells;
927 // maxAbsId - abs id of cell with highest energy
928 // e - cluster energy
930 // xr,yr,zr - x,y,z coordinates of cell with absId inside SM
932 // Shift index taking into account the difference between standard SM
933 // and SM of half size in phi direction
934 const Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
935 static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
936 static Int_t iphim, ietam;
937 static AliEMCALShishKebabTrd1Module *mod = 0;
940 static Int_t nSupModM, nModuleM, nIphiM, nIetaM, iphiM, ietaM;
941 static Int_t iphimM, ietamM, maxAbsIdCopy=-1;
942 static AliEMCALShishKebabTrd1Module *modM = 0;
943 static Double_t distCorr;
945 if(!CheckAbsCellId(absId)) return kFALSE;
947 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
948 GetModulePhiEtaIndexInSModule(nSupMod, nModule, iphim, ietam);
949 GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta);
950 mod = GetShishKebabModule(ietam);
952 if(absId != maxAbsId) {
954 if(maxAbsIdCopy != maxAbsId) {
955 GetCellIndex(maxAbsId, nSupModM, nModuleM, nIphiM, nIetaM);
956 GetModulePhiEtaIndexInSModule(nSupModM, nModuleM, iphimM, ietamM);
957 GetCellPhiEtaIndexInSModule(nSupModM,nModuleM,nIphiM,nIetaM, iphiM, ietaM);
958 modM = GetShishKebabModule(ietamM); // do I need this ?
959 maxAbsIdCopy = maxAbsId;
962 distCorr = GetEtaModuleSize()*(ietam-ietamM)/TMath::Tan(modM->GetTheta()); // Stay here
963 //printf(" distCorr %f | dist %f | ietam %i -> etamM %i\n", distCorr, dist, ietam, ietamM);
965 // distEff += distCorr;
967 // Bad resolution in this case, strong bias vs phi
969 mod->GetPositionAtCenterCellLine(nIeta, distEff, v); // Stay here
970 xr = v.Y() - fParSM[0];
971 zr = v.X() - fParSM[2];
974 yr = fCentersOfCellsPhiDir.At(iphi);
976 yr = fCentersOfCellsPhiDir.At(iphi + kphiIndexShift);
978 AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
983 //________________________________________________________________________________________________
984 void AliEMCALGeometry::CreateListOfTrd1Modules()
986 // Generate the list of Trd1 modules
987 // which will make up the EMCAL
990 AliDebug(2,Form(" AliEMCALGeometry::CreateListOfTrd1Modules() started "));
992 AliEMCALShishKebabTrd1Module *mod=0, *mTmp=0; // current module
993 if(fShishKebabTrd1Modules == 0) {
994 fShishKebabTrd1Modules = new TList;
995 fShishKebabTrd1Modules->SetName("ListOfTRD1");
996 for(int iz=0; iz< GetNZ(); iz++) {
998 mod = new AliEMCALShishKebabTrd1Module(TMath::Pi()/2.,this);
1000 mTmp = new AliEMCALShishKebabTrd1Module(*mod);
1003 fShishKebabTrd1Modules->Add(mod);
1006 AliDebug(2,Form(" Already exits : "));
1008 mod = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(fShishKebabTrd1Modules->GetSize()-1);
1009 fEtaMaxOfTRD1 = mod->GetMaxEtaOfModule(0);
1011 AliDebug(2,Form(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n",
1012 fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1));
1014 // Jun 01, 2006 - ALICE numbering scheme
1015 // define grid for cells in eta(z) and x directions in local coordinates system of SM
1016 // Works just for 2x2 case only -- ?? start here
1019 // Define grid for cells in phi(y) direction in local coordinates system of SM
1020 // as for 2X2 as for 3X3 - Nov 8,2006
1022 AliDebug(2,Form(" Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize()));
1023 Int_t ind=0; // this is phi index
1024 Int_t ieta=0, nModule=0, iphiTemp;
1025 Double_t xr, zr, theta, phi, eta, r, x,y;
1027 Double_t ytCenterModule=0.0, ytCenterCell=0.0;
1029 fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
1030 fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
1032 Double_t r0 = GetIPDistance() + GetLongModuleSize()/2.;
1033 for(Int_t it=0; it<fNPhi; it++) { // cycle on modules
1034 ytCenterModule = -fParSM[1] + fPhiModuleSize*(2*it+1)/2; // center of module
1035 for(Int_t ic=0; ic<fNPHIdiv; ic++) { // cycle on cells in module
1037 ytCenterCell = ytCenterModule + fPhiTileSize *(2*ic-1)/2.;
1038 } else if(fNPHIdiv==3){
1039 ytCenterCell = ytCenterModule + fPhiTileSize *(ic-1);
1040 } else if(fNPHIdiv==1){
1041 ytCenterCell = ytCenterModule;
1043 fCentersOfCellsPhiDir.AddAt(ytCenterCell,ind);
1044 // Define grid on phi direction
1045 // Grid is not the same for different eta bin;
1046 // Effect is small but is still here
1047 phi = TMath::ATan2(ytCenterCell, r0);
1048 fPhiCentersOfCells.AddAt(phi, ind);
1050 AliDebug(2,Form(" ind %2.2i : y %8.3f ", ind, fCentersOfCellsPhiDir.At(ind)));
1055 fCentersOfCellsEtaDir.Set(fNZ *fNETAdiv);
1056 fCentersOfCellsXDir.Set(fNZ *fNETAdiv);
1057 fEtaCentersOfCells.Set(fNZ *fNETAdiv * fNPhi*fNPHIdiv);
1058 AliDebug(2,Form(" Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize()));
1059 for(Int_t it=0; it<fNZ; it++) {
1060 AliEMCALShishKebabTrd1Module *trd1 = GetShishKebabModule(it);
1062 for(Int_t ic=0; ic<fNETAdiv; ic++) {
1064 trd1->GetCenterOfCellInLocalCoordinateofSM(ic, xr, zr); // case of 2X2
1065 GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1067 trd1->GetCenterOfCellInLocalCoordinateofSM_3X3(ic, xr, zr); // case of 3X3
1068 GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1070 trd1->GetCenterOfCellInLocalCoordinateofSM_1X1(xr, zr); // case of 1X1
1071 GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta);
1073 fCentersOfCellsXDir.AddAt(float(xr) - fParSM[0],ieta);
1074 fCentersOfCellsEtaDir.AddAt(float(zr) - fParSM[2],ieta);
1075 // Define grid on eta direction for each bin in phi
1076 for(int iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
1077 x = xr + trd1->GetRadius();
1078 y = fCentersOfCellsPhiDir[iphi];
1079 r = TMath::Sqrt(x*x + y*y + zr*zr);
1080 theta = TMath::ACos(zr/r);
1081 eta = AliEMCALShishKebabTrd1Module::ThetaToEta(theta);
1082 // ind = ieta*fCentersOfCellsPhiDir.GetSize() + iphi;
1083 ind = iphi*fCentersOfCellsEtaDir.GetSize() + ieta;
1084 fEtaCentersOfCells.AddAt(eta, ind);
1086 //printf(" ieta %i : xr + trd1->GetRadius() %f : zr %f : eta %f \n", ieta, xr + trd1->GetRadius(), zr, eta);
1089 for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1090 AliDebug(2,Form(" ind %2.2i : z %8.3f : x %8.3f", i+1,
1091 fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i)));
1096 //________________________________________________________________________________________________
1097 void AliEMCALGeometry::GetGlobal(const Double_t *loc, Double_t *glob, int ind) const
1099 // Figure out the global numbering
1100 // of a given supermodule from the
1101 // local numbering and the transformation
1102 // matrix stored by the geometry manager (allows for misaligned
1105 if(ind>=0 && ind < GetNumberOfSuperModules()) {
1106 TString volpath = "ALIC_1/XEN1_1/SMOD_";
1109 if(GetKey110DEG() && ind>=10) {
1110 volpath = "ALIC_1/XEN1_1/SM10_";
1111 volpath += ind-10+1;
1114 if(!gGeoManager->cd(volpath.Data()))
1115 AliFatal(Form("AliEMCALGeometry::GeoManager cannot find path %s!",volpath.Data()));
1117 TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1119 m->LocalToMaster(loc, glob);
1121 AliFatal("Geo matrixes are not loaded \n") ;
1126 //________________________________________________________________________________________________
1127 void AliEMCALGeometry::GetGlobal(const TVector3 &vloc, TVector3 &vglob, int ind) const
1129 //Figure out the global numbering
1130 //of a given supermodule from the
1131 //local numbering given a 3-vector location
1133 static Double_t tglob[3], tloc[3];
1135 GetGlobal(tloc, tglob, ind);
1136 vglob.SetXYZ(tglob[0], tglob[1], tglob[2]);
1139 //________________________________________________________________________________________________
1140 void AliEMCALGeometry::GetGlobal(Int_t absId , double glob[3]) const
1142 // Alice numbering scheme - Jun 03, 2006
1143 static Int_t nSupMod, nModule, nIphi, nIeta;
1144 static double loc[3];
1146 if (!gGeoManager || !gGeoManager->IsClosed()) {
1147 AliError("Can't get the global coordinates! gGeoManager doesn't exist or it is still open!");
1151 glob[0]=glob[1]=glob[2]=0.0; // bad case
1152 if(RelPosCellInSModule(absId, loc)) {
1153 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1155 TString volpath = "ALIC_1/XEN1_1/SMOD_";
1156 volpath += (nSupMod+1);
1158 if(GetKey110DEG() && nSupMod>=10) {
1159 volpath = "ALIC_1/XEN1_1/SM10_";
1160 volpath += (nSupMod-10+1);
1162 if(!gGeoManager->cd(volpath.Data()))
1163 AliFatal(Form("GeoManager cannot find path %s!",volpath.Data()));
1165 TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1167 m->LocalToMaster(loc, glob);
1169 AliFatal("Geo matrixes are not loaded \n") ;
1174 //___________________________________________________________________
1175 void AliEMCALGeometry::GetGlobal(Int_t absId , TVector3 &vglob) const
1177 // Alice numbering scheme - Jun 03, 2006
1178 static Double_t glob[3];
1180 GetGlobal(absId, glob);
1181 vglob.SetXYZ(glob[0], glob[1], glob[2]);
1185 //____________________________________________________________________________
1186 void AliEMCALGeometry::GetGlobal(const AliRecPoint* /*rp*/, TVector3& /* vglob */) const
1188 AliFatal(Form("Please use GetGlobalEMCAL(recPoint,gpos) instead of GetGlobal!"));
1191 //_________________________________________________________________________________
1192 void AliEMCALGeometry::GetGlobalEMCAL(const AliEMCALRecPoint *rp, TVector3 &vglob) const
1194 // Figure out the global numbering
1195 // of a given supermodule from the
1196 // local numbering for RecPoints
1198 static TVector3 vloc;
1199 static Int_t nSupMod, nModule, nIphi, nIeta;
1201 const AliEMCALRecPoint *rpTmp = rp;
1202 const AliEMCALRecPoint *rpEmc = rpTmp;
1204 GetCellIndex(rpEmc->GetAbsId(0), nSupMod, nModule, nIphi, nIeta);
1205 rpTmp->GetLocalPosition(vloc);
1206 GetGlobal(vloc, vglob, nSupMod);
1209 //________________________________________________________________________________________________
1210 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Double_t &eta,Double_t &phi) const
1212 // Nov 16, 2006- float to double
1213 // version for TRD1 only
1214 static TVector3 vglob;
1215 GetGlobal(absId, vglob);
1220 //________________________________________________________________________________________________
1221 void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Float_t &eta,Float_t &phi) const
1223 // Nov 16,2006 - should be discard in future
1224 static TVector3 vglob;
1225 GetGlobal(absId, vglob);
1226 eta = float(vglob.Eta());
1227 phi = float(vglob.Phi());
1230 //________________________________________________________________________________________________
1231 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSM(Int_t nSupMod, Double_t &phiMin, Double_t &phiMax) const
1233 // 0<= nSupMod <=11; phi in rad
1235 if(nSupMod<0 || nSupMod >11) return kFALSE;
1237 phiMin = fPhiBoundariesOfSM[2*i];
1238 phiMax = fPhiBoundariesOfSM[2*i+1];
1242 //________________________________________________________________________________________________
1243 Bool_t AliEMCALGeometry::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const
1245 // 0<= nPhiSec <=4; phi in rad
1246 // 0; gap boundaries between 0th&2th | 1th&3th SM
1247 // 1; gap boundaries between 2th&4th | 3th&5th SM
1248 // 2; gap boundaries between 4th&6th | 5th&7th SM
1249 // 3; gap boundaries between 6th&8th | 7th&9th SM
1250 // 4; gap boundaries between 8th&10th | 9th&11th SM
1251 if(nPhiSec<0 || nPhiSec >4) return kFALSE;
1252 phiMin = fPhiBoundariesOfSM[2*nPhiSec+1];
1253 phiMax = fPhiBoundariesOfSM[2*nPhiSec+2];
1257 //________________________________________________________________________________________________
1258 Bool_t AliEMCALGeometry::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const
1260 // Return false if phi belongs a phi cracks between SM
1264 if(TMath::Abs(eta) > fEtaMaxOfTRD1) return kFALSE;
1266 phi = TVector2::Phi_0_2pi(phi); // move phi to (0,2pi) boundaries
1267 for(i=0; i<6; i++) {
1268 if(phi>=fPhiBoundariesOfSM[2*i] && phi<=fPhiBoundariesOfSM[2*i+1]) {
1270 if(eta < 0.0) nSupMod++;
1271 AliDebug(1,Form("eta %f phi %f(%5.2f) : nSupMod %i : #bound %i", eta,phi,phi*TMath::RadToDeg(), nSupMod,i));
1278 //________________________________________________________________________________________________
1279 Bool_t AliEMCALGeometry::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi, Int_t &absId) const
1282 // stay here - phi problem as usual
1283 static Int_t nSupMod, i, ieta, iphi, etaShift, nphi;
1284 static Double_t absEta=0.0, d=0.0, dmin=0.0, phiLoc;
1285 absId = nSupMod = - 1;
1286 if(SuperModuleNumberFromEtaPhi(eta, phi, nSupMod)) {
1288 phi = TVector2::Phi_0_2pi(phi);
1289 phiLoc = phi - fPhiCentersOfSM[nSupMod/2];
1290 nphi = fPhiCentersOfCells.GetSize();
1292 phiLoc = phi - 190.*TMath::DegToRad();
1296 dmin = TMath::Abs(fPhiCentersOfCells[0]-phiLoc);
1298 for(i=1; i<nphi; i++) {
1299 d = TMath::Abs(fPhiCentersOfCells[i] - phiLoc);
1304 // printf(" i %i : d %f : dmin %f : fPhiCentersOfCells[i] %f \n", i, d, dmin, fPhiCentersOfCells[i]);
1306 // odd SM are turned with respect of even SM - reverse indexes
1307 AliDebug(2,Form(" iphi %i : dmin %f (phi %f, phiLoc %f ) ", iphi, dmin, phi, phiLoc));
1309 absEta = TMath::Abs(eta);
1310 etaShift = iphi*fCentersOfCellsEtaDir.GetSize();
1311 dmin = TMath::Abs(fEtaCentersOfCells[etaShift]-absEta);
1313 for(i=1; i<fCentersOfCellsEtaDir.GetSize(); i++) {
1314 d = TMath::Abs(fEtaCentersOfCells[i+etaShift] - absEta);
1320 AliDebug(2,Form(" ieta %i : dmin %f (eta=%f) : nSupMod %i ", ieta, dmin, eta, nSupMod));
1322 if(eta<0) iphi = (nphi-1) - iphi;
1323 absId = GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);
1330 //________________________________________________________________________________________________
1331 AliEMCALShishKebabTrd1Module* AliEMCALGeometry::GetShishKebabModule(Int_t neta) const
1333 //This method was too long to be
1334 //included in the header file - the
1335 //rule checker complained about it's
1336 //length, so we move it here. It returns the
1337 //shishkebabmodule at a given eta index point.
1339 static AliEMCALShishKebabTrd1Module* trd1=0;
1340 if(fShishKebabTrd1Modules && neta>=0 && neta<fShishKebabTrd1Modules->GetSize()) {
1341 trd1 = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(neta);
1346 //________________________________________________________________________________________________
1347 Int_t AliEMCALGeometry::GetAbsTRUNumberFromNumberInSm(const Int_t row, const Int_t col, const Int_t sm)
1349 Int_t itru = row + col*GetNModulesInTRUPhi() + sm*GetNTRU();
1350 // printf(" GetAbsTRUNumberFromNumberInSm : row %2i col %2i sm %2i -> itru %2i\n", row, col, sm, itru);
1354 //________________________________________________________________________________________________
1355 void AliEMCALGeometry::Browse(TBrowser* b)
1357 //Browse the modules
1358 if(fShishKebabTrd1Modules) b->Add(fShishKebabTrd1Modules);
1361 //________________________________________________________________________________________________
1362 Bool_t AliEMCALGeometry::IsFolder() const
1364 //Check if fShishKebabTrd1Modules is in folder
1365 if(fShishKebabTrd1Modules) return kTRUE;
1369 //________________________________________________________________________________________________
1370 Double_t AliEMCALGeometry::GetPhiCenterOfSM(Int_t nsupmod) const
1372 //returns center of supermodule in phi
1374 return fPhiCentersOfSM[i];
1377 //____________________________________________________________________________
1378 Bool_t AliEMCALGeometry::Impact(const TParticle * particle) const
1380 // Tells if a particle enters EMCAL
1383 TVector3 vtx(particle->Vx(),particle->Vy(),particle->Vz());
1384 TVector3 vimpact(0,0,0);
1385 ImpactOnEmcal(vtx,particle->Theta(),particle->Phi(),AbsID,vimpact);
1390 //____________________________________________________________________________
1391 void AliEMCALGeometry::ImpactOnEmcal(TVector3 vtx, Double_t theta, Double_t phi,
1392 Int_t & absId, TVector3 & vimpact) const
1394 // calculates the impact coordinates on EMCAL (centre of a tower/not on EMCAL surface)
1395 // of a neutral particle
1396 // emitted in the vertex vtx[3] with direction theta and phi in the ALICE global coordinate system
1398 TVector3 p(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) ;
1400 vimpact.SetXYZ(0,0,0);
1402 if(phi==0 || theta==0) return;
1405 Double_t factor = (GetIPDistance()-vtx[1])/p[1];
1406 direction = vtx + factor*p;
1409 AliFatal("Geo manager not initialized\n");
1411 //from particle direction -> tower hitted
1412 GetAbsCellIdFromEtaPhi(direction.Eta(),direction.Phi(),absId);
1414 //tower absID hitted -> tower/module plane (evaluated at the center of the tower)
1415 Int_t nSupMod, nModule, nIphi, nIeta;
1416 Double_t loc[3],loc2[3],loc3[3];
1417 Double_t glob[3]={},glob2[3]={},glob3[3]={};
1419 if(!RelPosCellInSModule(absId,loc)) return;
1421 //loc is cell center of tower
1422 GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1424 //look at 2 neighbours-s cell using nIphi={0,1} and nIeta={0,1}
1425 Int_t nIphi2,nIeta2,absId2,absId3;
1426 if(nIeta==0) nIeta2=1;
1428 absId2=GetAbsCellId(nSupMod,nModule,nIphi,nIeta2);
1429 if(nIphi==0) nIphi2=1;
1431 absId3=GetAbsCellId(nSupMod,nModule,nIphi2,nIeta);
1433 //2nd point on emcal cell plane
1434 if(!RelPosCellInSModule(absId2,loc2)) return;
1436 //3rd point on emcal cell plane
1437 if(!RelPosCellInSModule(absId3,loc3)) return;
1439 TString volpath = "ALIC_1/XEN1_1/SMOD_";
1440 volpath += (nSupMod+1);
1442 if(GetKey110DEG() && nSupMod>=10) {
1443 volpath = "ALIC_1/XEN1_1/SM10_";
1444 volpath += (nSupMod-10+1);
1446 if(!gGeoManager->cd(volpath.Data())){
1447 AliFatal(Form("GeoManager cannot find path %s!",volpath.Data()))
1450 TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
1452 m->LocalToMaster(loc, glob);
1453 m->LocalToMaster(loc2, glob2);
1454 m->LocalToMaster(loc3, glob3);
1456 AliFatal("Geo matrixes are not loaded \n") ;
1459 //Equation of Plane from glob,glob2,glob3 (Ax+By+Cz+D=0)
1460 Double_t A = glob[1]*(glob2[2]-glob3[2]) + glob2[1]*(glob3[2]-glob[2]) + glob3[1]*(glob[2]-glob2[2]);
1461 Double_t B = glob[2]*(glob2[0]-glob3[0]) + glob2[2]*(glob3[0]-glob[0]) + glob3[2]*(glob[0]-glob2[0]);
1462 Double_t C = glob[0]*(glob2[1]-glob3[1]) + glob2[0]*(glob3[1]-glob[1]) + glob3[0]*(glob[1]-glob2[1]);
1463 Double_t D = glob[0]*(glob2[1]*glob3[2]-glob3[1]*glob2[2]) + glob2[0]*(glob3[1]*glob[2]-glob[1]*glob3[2]) + glob3[0]*(glob[1]*glob2[2]-glob2[1]*glob[2]);
1466 //shift equation of plane from tower/module center to surface along vector (A,B,C) normal to tower/module plane
1467 Double_t dist = GetLongModuleSize()/2.;
1468 Double_t norm = TMath::Sqrt(A*A+B*B+C*C);
1469 Double_t glob4[3]={};
1470 TVector3 dir(A,B,C);
1471 TVector3 point(glob[0],glob[1],glob[2]);
1472 if(point.Dot(dir)<0) dist*=-1;
1473 glob4[0]=glob[0]-dist*A/norm;
1474 glob4[1]=glob[1]-dist*B/norm;
1475 glob4[2]=glob[2]-dist*C/norm;
1476 D = glob4[0]*A + glob4[1]*B + glob4[2]*C ;
1479 //Line determination (2 points for equation of line : vtx and direction)
1480 //impact between line (particle) and plane (module/tower plane)
1481 Double_t den = A*(vtx(0)-direction(0)) + B*(vtx(1)-direction(1)) + C*(vtx(2)-direction(2));
1483 printf("ImpactOnEmcal() No solution :\n");
1487 Double_t length = A*vtx(0)+B*vtx(1)+C*vtx(2)+D;
1490 vimpact.SetXYZ(vtx(0)+length*(direction(0)-vtx(0)),vtx(1)+length*(direction(1)-vtx(1)),vtx(2)+length*(direction(2)-vtx(2)));
1492 //shift vimpact from tower/module surface to center along vector (A,B,C) normal to tower/module plane
1493 vimpact.SetXYZ(vimpact(0)+dist*A/norm,vimpact(1)+dist*B/norm,vimpact(2)+dist*C/norm);