new function EtaPhiFromIndex
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALGeometry.cxx
CommitLineData
e5a93224 1/**************************************************************************
2012850d 2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$*/
17
18//_________________________________________________________________________
19// Geometry class for EMCAL : singleton
b13bbe81 20// EMCAL consists of layers of scintillator and lead
ffa6d63b 21// Places the the Barrel Geometry of The EMCAL at Midrapidity
d87bd045 22// between 80 and 180(or 190) degrees of Phi and
ffa6d63b 23// -0.7 to 0.7 in eta
24// Number of Modules and Layers may be controlled by
25// the name of the instance defined
b13bbe81 26//*-- Author: Sahal Yacoob (LBL / UCT)
27// and : Yves Schutz (SUBATECH)
28// and : Jennifer Klay (LBL)
1963b290 29// SHASHLYK : Aleksei Pavlinov (WSU)
c63c3c5d 30// SuperModules -> module(or tower) -> cell
2012850d 31
2012850d 32// --- AliRoot header files ---
1ceeec56 33#include <assert.h>
e52475ed 34#include "Riostream.h"
35
ca8f5bd0 36#include <TMath.h>
116cbefd 37#include <TVector3.h>
dc7da436 38 //#include <TArrayD.h>
c63c3c5d 39#include <TObjArray.h>
e52475ed 40#include <TGeoManager.h>
41#include <TGeoNode.h>
42#include <TGeoMatrix.h>
f0377b23 43#include <TMatrixD.h>
d434833b 44#include <TObjString.h>
f0377b23 45#include <TClonesArray.h>
173558f2 46
ca8f5bd0 47// -- ALICE Headers.
d64c959b 48//#include "AliConst.h"
e5a93224 49#include "AliLog.h"
173558f2 50
ca8f5bd0 51// --- EMCAL headers
52#include "AliEMCALGeometry.h"
e52475ed 53#include "AliEMCALShishKebabTrd1Module.h"
e52475ed 54#include "AliEMCALRecPoint.h"
f0377b23 55#include "AliEMCALDigit.h"
d434833b 56#include "AliEMCALHistoUtilities.h"
57#include "AliEMCALAlignData.h"
2012850d 58
925e6570 59ClassImp(AliEMCALGeometry)
2012850d 60
d434833b 61// these initialisations are needed for a singleton
62AliEMCALGeometry *AliEMCALGeometry::fgGeom = 0;
63Bool_t AliEMCALGeometry::fgInit = kFALSE;
64AliEMCALAlignData *AliEMCALGeometry::fgAlignData = 0;
65
dc7da436 66
67
68AliEMCALGeometry::AliEMCALGeometry() : AliGeometry()
69{
70 // default ctor only for internal usage (singleton)
71 // must be kept public for root persistency purposes, but should never be called by the outside world
72 // CreateListOfTrd1Modules();
73 AliDebug(2, "AliEMCALGeometry : default ctor ");
74}
75//______________________________________________________________________
76AliEMCALGeometry::AliEMCALGeometry(const Text_t* name, const Text_t* title) :
77AliGeometry(name, title) {// ctor only for internal usage (singleton)
78 AliDebug(2, Form("AliEMCALGeometry(%s,%s) ", name,title));
79 Init();
80 CreateListOfTrd1Modules();
81}
b13bbe81 82//______________________________________________________________________
83AliEMCALGeometry::~AliEMCALGeometry(void){
84 // dtor
2012850d 85}
395c7ba2 86//______________________________________________________________________
87void AliEMCALGeometry::Init(void){
88 // Initializes the EMCAL parameters
fdebddeb 89 // naming convention : GUV_WX_N_ gives the composition of a tower
395c7ba2 90 // WX inform about the composition of the EM calorimeter section:
fdebddeb 91 // thickness in mm of Pb radiator (W) and of scintillator (X), and number of scintillator layers (N)
92 // New geometry: EMCAL_55_25
1963b290 93 // 24-aug-04 for shish-kebab
94 // SHISH_25 or SHISH_62
c63c3c5d 95 // 11-oct-05 - correction for pre final design
96 // Feb 06,2006 - decrease the weight of EMCAL
fc575e27 97
b44d5aa4 98 fAdditionalOpts[0] = "nl="; // number of sampling layers (fNECLayers)
99 fAdditionalOpts[1] = "pbTh="; // cm, Thickness of the Pb (fECPbRadThick)
100 fAdditionalOpts[2] = "scTh="; // cm, Thickness of the Sc (fECScintThick)
101 fAdditionalOpts[3] = "latSS="; // cm, Thickness of lateral steel strip (fLateralSteelStrip)
fc575e27 102
103 fNAdditionalOpts = sizeof(fAdditionalOpts) / sizeof(char*);
104
fdebddeb 105 fgInit = kFALSE; // Assume failed until proven otherwise.
fc575e27 106 fGeoName = GetName();
107 fGeoName.ToUpper();
d87bd045 108 fKey110DEG = 0;
fc575e27 109 if(fGeoName.Contains("110DEG")) fKey110DEG = 1; // for GetAbsCellId
e52475ed 110 fShishKebabTrd1Modules = 0;
111 fTrd2AngleY = f2Trd2Dy2 = fEmptySpace = fTubsR = fTubsTurnAngle = 0;
1963b290 112
113 fNZ = 114; // granularity along Z (eta)
114 fNPhi = 168; // granularity in phi (azimuth)
115 fArm1PhiMin = 60.0; // degrees, Starting EMCAL Phi position
116 fArm1PhiMax = 180.0; // degrees, Ending EMCAL Phi position
117 fArm1EtaMin = -0.7; // pseudorapidity, Starting EMCAL Eta position
118 fArm1EtaMax = +0.7; // pseudorapidity, Ending EMCAL Eta position
119 fIPDistance = 454.0; // cm, Radial distance to inner surface of EMCAL
905263da 120 fPhiGapForSM = 0.; // cm, only for final TRD1 geometry
e52475ed 121 for(int i=0; i<12; i++) fMatrixOfSM[i] = 0;
1963b290 122
123 // geometry
fc575e27 124 if(fGeoName.Contains("SHISH")){ // Only shahslyk now
905263da 125 // 7-sep-05; integration issue
126 fArm1PhiMin = 80.0; // 60 -> 80
c01485dd 127 fArm1PhiMax = 180.0; // 180 -> 190
905263da 128
129 fNumberOfSuperModules = 10; // 12 = 6 * 2 (6 in phi, 2 in Z);
1963b290 130 fSteelFrontThick = 2.54; // 9-sep-04
131 fIPDistance = 460.0;
132 fFrontSteelStrip = fPassiveScintThick = 0.0; // 13-may-05
133 fLateralSteelStrip = 0.025; // before MAY 2005
134 fPhiModuleSize = fEtaModuleSize = 11.4;
135 fPhiTileSize = fEtaTileSize = 5.52; // (11.4-5.52*2)/2. = 0.18 cm (wall thickness)
136 fNPhi = 14;
137 fNZ = 30;
138 fAlFrontThick = fGap2Active = 0;
139 fNPHIdiv = fNETAdiv = 2;
140
141 fNECLayers = 62;
142 fECScintThick = fECPbRadThickness = 0.2;
143 fSampling = 1.; // 30-aug-04 - should be calculated
fc575e27 144 if(fGeoName.Contains("TWIST")) { // all about EMCAL module
1963b290 145 fNZ = 27; // 16-sep-04
fc575e27 146 } else if(fGeoName.Contains("TRD")) {
1963b290 147 fIPDistance = 428.0; // 11-may-05
148 fSteelFrontThick = 0.0; // 3.17 -> 0.0; 28-mar-05 : no stell plate
149 fNPhi = 12;
150 fSampling = 12.327;
151 fPhiModuleSize = fEtaModuleSize = 12.26;
152 fNZ = 26; // 11-oct-04
153 fTrd1Angle = 1.3; // in degree
154// 18-nov-04; 1./0.08112=12.327
155// http://pdsfweb01.nersc.gov/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
fc575e27 156 if(fGeoName.Contains("TRD1")) { // 30-jan-05
1963b290 157 // for final design
905263da 158 fPhiGapForSM = 2.; // cm, only for final TRD1 geometry
fc575e27 159 if(fGeoName.Contains("MAY05") || fGeoName.Contains("WSUC") || fGeoName.Contains("FINAL")){
1963b290 160 fNumberOfSuperModules = 12; // 20-may-05
fc575e27 161 if(fGeoName.Contains("WSUC")) fNumberOfSuperModules = 1; // 27-may-05
1963b290 162 fNECLayers = 77; // (13-may-05 from V.Petrov)
163 fPhiModuleSize = 12.5; // 20-may-05 - rectangular shape
164 fEtaModuleSize = 11.9;
165 fECScintThick = fECPbRadThickness = 0.16;// (13-may-05 from V.Petrov)
166 fFrontSteelStrip = 0.025;// 0.025cm = 0.25mm (13-may-05 from V.Petrov)
167 fLateralSteelStrip = 0.01; // 0.01cm = 0.1mm (13-may-05 from V.Petrov) - was 0.025
168 fPassiveScintThick = 0.8; // 0.8cm = 8mm (13-may-05 from V.Petrov)
169 fNZ = 24;
170 fTrd1Angle = 1.5; // 1.3 or 1.5
905263da 171
fc575e27 172 if(fGeoName.Contains("FINAL")) { // 9-sep-05
905263da 173 fNumberOfSuperModules = 10;
fc575e27 174 if(fGeoName.Contains("110DEG")) {
d87bd045 175 fNumberOfSuperModules = 12;// last two modules have size 10 degree in phi (180<phi<190)
176 fArm1PhiMax = 200.0; // for XEN1 and turn angle of super modules
177 }
905263da 178 fPhiModuleSize = 12.26 - fPhiGapForSM / Float_t(fNPhi); // first assumption
179 fEtaModuleSize = fPhiModuleSize;
fc575e27 180 if(fGeoName.Contains("HUGE")) fNECLayers *= 3; // 28-oct-05 for analysing leakage
905263da 181 }
1963b290 182 }
fc575e27 183 } else if(fGeoName.Contains("TRD2")) { // 30-jan-05
1963b290 184 fSteelFrontThick = 0.0; // 11-mar-05
185 fIPDistance+= fSteelFrontThick; // 1-feb-05 - compensate absence of steel plate
186 fTrd1Angle = 1.64; // 1.3->1.64
187 fTrd2AngleY = fTrd1Angle; // symmetric case now
188 fEmptySpace = 0.2; // 2 mm
189 fTubsR = fIPDistance; // 31-jan-05 - as for Fred case
190
191 fPhiModuleSize = fTubsR*2.*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.);
192 fPhiModuleSize -= fEmptySpace/2.; // 11-mar-05
193 fEtaModuleSize = fPhiModuleSize; // 20-may-05
194 fTubsTurnAngle = 3.;
195 }
196 fNPHIdiv = fNETAdiv = 2; // 13-oct-04 - division again
fc575e27 197 if(fGeoName.Contains("3X3")) { // 23-nov-04
1963b290 198 fNPHIdiv = fNETAdiv = 3;
fc575e27 199 } else if(fGeoName.Contains("4X4")) {
1963b290 200 fNPHIdiv = fNETAdiv = 4;
201 }
202 }
fc575e27 203 if(fGeoName.Contains("25")){
1963b290 204 fNECLayers = 25;
205 fECScintThick = fECPbRadThickness = 0.5;
206 }
fc575e27 207 if(fGeoName.Contains("WSUC")){ // 18-may-05 - about common structure
1963b290 208 fShellThickness = 30.; // should be change
209 fNPhi = fNZ = 4;
210 }
c63c3c5d 211
fc575e27 212 CheckAdditionalOptions();
c63c3c5d 213
31b39a2e 214 fPhiTileSize = fPhiModuleSize/2. - fLateralSteelStrip; // 13-may-05
215 fEtaTileSize = fEtaModuleSize/2. - fLateralSteelStrip; // 13-may-05
216
1963b290 217 // constant for transition absid <--> indexes
218 fNCellsInTower = fNPHIdiv*fNETAdiv;
219 fNCellsInSupMod = fNCellsInTower*fNPhi*fNZ;
220 fNCells = fNCellsInSupMod*fNumberOfSuperModules;
fc575e27 221 if(fGeoName.Contains("110DEG")) fNCells -= fNCellsInSupMod;
1963b290 222
223 fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness);
fc575e27 224 if(fGeoName.Contains("MAY05")) fLongModuleSize += (fFrontSteelStrip + fPassiveScintThick);
1963b290 225
226 // 30-sep-04
fc575e27 227 if(fGeoName.Contains("TRD")) {
1963b290 228 f2Trd1Dx2 = fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.);
fc575e27 229 if(fGeoName.Contains("TRD2")) { // 27-jan-05
1963b290 230 f2Trd2Dy2 = fPhiModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.);
231 }
232 }
fc575e27 233 } else Fatal("Init", "%s is an undefined geometry!", fGeoName.Data()) ;
fdebddeb 234
1963b290 235 fNPhiSuperModule = fNumberOfSuperModules/2;
236 if(fNPhiSuperModule<1) fNPhiSuperModule = 1;
fdebddeb 237 //There is always one more scintillator than radiator layer because of the first block of aluminium
238 fShellThickness = fAlFrontThick + fGap2Active + fNECLayers*GetECScintThick()+(fNECLayers-1)*GetECPbRadThick();
fc575e27 239 if(fGeoName.Contains("SHISH")) {
1963b290 240 fShellThickness = fSteelFrontThick + fLongModuleSize;
fc575e27 241 if(fGeoName.Contains("TWIST")) { // 13-sep-04
1963b290 242 fShellThickness = TMath::Sqrt(fLongModuleSize*fLongModuleSize + fPhiModuleSize*fEtaModuleSize);
243 fShellThickness += fSteelFrontThick;
fc575e27 244 } else if(fGeoName.Contains("TRD")) { // 1-oct-04
1963b290 245 fShellThickness = TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2);
246 fShellThickness += fSteelFrontThick;
e52475ed 247 // Local coordinates
248 fParSM[0] = GetShellThickness()/2.;
249 fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
250 fParSM[2] = 350./2.;
1963b290 251 }
252 }
fdebddeb 253
395c7ba2 254 fZLength = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax); // Z coverage
255 fEnvelop[0] = fIPDistance; // mother volume inner radius
256 fEnvelop[1] = fIPDistance + fShellThickness; // mother volume outer r.
257 fEnvelop[2] = 1.00001*fZLength; // add some padding for mother volume.
258
d434833b 259 if(fgAlignData != NULL) {
260 // Number of modules is read from Alignment DB if exists
261 fNumberOfSuperModules = fgAlignData->GetNSuperModules();
262 }
263
395c7ba2 264 fgInit = kTRUE;
265
e5a93224 266 if (AliDebugLevel()>=2) {
fc575e27 267 printf("Init: geometry of EMCAL named %s is as follows:\n", fGeoName.Data());
e52475ed 268 printf( " ECAL : %d x (%f cm Pb, %f cm Sc) \n",
269 GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ;
270 printf(" fSampling %5.2f \n", fSampling );
fc575e27 271 if(fGeoName.Contains("SHISH")){
1963b290 272 printf(" fIPDistance %6.3f cm \n", fIPDistance);
273 if(fSteelFrontThick>0.)
274 printf(" fSteelFrontThick %6.3f cm \n", fSteelFrontThick);
275 printf(" fNPhi %i | fNZ %i \n", fNPhi, fNZ);
d87bd045 276 printf(" fNCellsInTower %i : fNCellsInSupMod %i : fNCells %i\n",fNCellsInTower, fNCellsInSupMod, fNCells);
fc575e27 277 if(fGeoName.Contains("MAY05")){
1963b290 278 printf(" fFrontSteelStrip %6.4f cm (thickness of front steel strip)\n",
279 fFrontSteelStrip);
280 printf(" fLateralSteelStrip %6.4f cm (thickness of lateral steel strip)\n",
281 fLateralSteelStrip);
282 printf(" fPassiveScintThick %6.4f cm (thickness of front passive Sc tile)\n",
283 fPassiveScintThick);
284 }
c63c3c5d 285 printf(" X:Y module size %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize);
286 printf(" X:Y tile size %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize);
287 printf(" #of sampling layers %i(fNECLayers) \n", fNECLayers);
288 printf(" fLongModuleSize %6.3f cm \n", fLongModuleSize);
1963b290 289 printf(" #supermodule in phi direction %i \n", fNPhiSuperModule );
290 }
fc575e27 291 if(fGeoName.Contains("TRD")) {
1963b290 292 printf(" fTrd1Angle %7.4f\n", fTrd1Angle);
293 printf(" f2Trd1Dx2 %7.4f\n", f2Trd1Dx2);
fc575e27 294 if(fGeoName.Contains("TRD2")) {
1963b290 295 printf(" fTrd2AngleY %7.4f\n", fTrd2AngleY);
296 printf(" f2Trd2Dy2 %7.4f\n", f2Trd2Dy2);
905263da 297 printf(" fTubsR %7.2f cm\n", fTubsR);
1963b290 298 printf(" fTubsTurnAngle %7.4f\n", fTubsTurnAngle);
905263da 299 printf(" fEmptySpace %7.4f cm\n", fEmptySpace);
fc575e27 300 } else if(fGeoName.Contains("TRD1") && fGeoName.Contains("FINAL")){
e52475ed 301 printf("SM dimensions(TRD1) : dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n",
302 fParSM[0],fParSM[1],fParSM[2]);
905263da 303 printf(" fPhiGapForSM %7.4f cm \n", fPhiGapForSM);
fc575e27 304 if(fGeoName.Contains("110DEG"))printf(" Last two modules have size 10 degree in phi (180<phi<190)\n");
1963b290 305 }
306 }
88cb7938 307 printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
1963b290 308 printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f\n",
309 GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() );
88cb7938 310 }
f0377b23 311 //TRU parameters. These parameters values are not the final ones.
312 fNTRU = 3 ;
313 fNTRUEta = 3 ;
314 fNTRUPhi = 1 ;
2012850d 315}
173558f2 316
c63c3c5d 317//______________________________________________________________________
318
fc575e27 319void AliEMCALGeometry::CheckAdditionalOptions()
320{
321 // Feb 06,2006
322 //Additional options that
323 //can be used to select
324 //the specific geometry of
325 //EMCAL to run
326
c63c3c5d 327 fArrayOpts = new TObjArray;
fc575e27 328 Int_t nopt = AliEMCALHistoUtilities::ParseString(fGeoName, *fArrayOpts);
c63c3c5d 329 if(nopt==1) { // no aditional option(s)
330 fArrayOpts->Delete();
331 delete fArrayOpts;
332 fArrayOpts = 0;
333 return;
334 }
335 for(Int_t i=1; i<nopt; i++){
336 TObjString *o = (TObjString*)fArrayOpts->At(i);
337
338 TString addOpt = o->String();
339 Int_t indj=-1;
fc575e27 340 for(Int_t j=0; j<fNAdditionalOpts; j++) {
341 TString opt = fAdditionalOpts[j];
c63c3c5d 342 if(addOpt.Contains(opt,TString::kIgnoreCase)) {
343 indj = j;
344 break;
345 }
346 }
347 if(indj<0) {
e5a93224 348 AliDebug(2,Form("<E> option |%s| unavailable : ** look to the file AliEMCALGeometry.h **\n",
349 addOpt.Data()));
c63c3c5d 350 assert(0);
351 } else {
e5a93224 352 AliDebug(2,Form("<I> option |%s| is valid : number %i : |%s|\n",
353 addOpt.Data(), indj, fAdditionalOpts[indj]));
c63c3c5d 354 if (addOpt.Contains("NL=",TString::kIgnoreCase)) {// number of sampling layers
355 sscanf(addOpt.Data(),"NL=%i", &fNECLayers);
e5a93224 356 AliDebug(2,Form(" fNECLayers %i (new) \n", fNECLayers));
b44d5aa4 357 } else if(addOpt.Contains("PBTH=",TString::kIgnoreCase)) {//Thickness of the Pb(fECPbRadThicknes)
c63c3c5d 358 sscanf(addOpt.Data(),"PBTH=%f", &fECPbRadThickness);
b44d5aa4 359 } else if(addOpt.Contains("SCTH=",TString::kIgnoreCase)) {//Thickness of the Sc(fECScintThick)
c63c3c5d 360 sscanf(addOpt.Data(),"SCTH=%f", &fECScintThick);
b44d5aa4 361 } else if(addOpt.Contains("LATSS=",TString::kIgnoreCase)) {// Thickness of lateral steel strip (fLateralSteelStrip)
362 sscanf(addOpt.Data(),"LATSS=%f", &fLateralSteelStrip);
363 AliDebug(2,Form(" fLateralSteelStrip %f (new) \n", fLateralSteelStrip));
c63c3c5d 364 }
365 }
366 }
367}
368
f0377b23 369//____________________________________________________________________________
356fd0a9 370void AliEMCALGeometry::FillTRU(const TClonesArray * digits, TClonesArray * ampmatrix, TClonesArray * timeRmatrix) {
f0377b23 371
372
356fd0a9 373// Orders digits ampitudes list in fNTRU TRUs (384 cells) per supermodule.
374// Each TRU is a TMatrixD, and they are kept in TClonesArrays. The number of
375// TRU in phi is fNTRUPhi, and the number of TRU in eta is fNTRUEta.
376// Last 2 modules are half size in Phi, I considered that the number of TRU
377// is maintained for the last modules but decision not taken. If different,
378// then this must be changed.
379
f0377b23 380
381 //Check data members
382
383 if(fNTRUEta*fNTRUPhi != fNTRU)
384 Error("FillTRU"," Wrong number of TRUS per Eta or Phi");
385
356fd0a9 386 //Initilize and declare variables
f0377b23 387 //List of TRU matrices initialized to 0.
356fd0a9 388 Int_t nCellsPhi = fNPhi*2/fNTRUPhi;
389 Int_t nCellsPhi2 = fNPhi/fNTRUPhi; //HalfSize modules
390 Int_t nCellsEta = fNZ*2/fNTRUEta;
f0377b23 391 Int_t id = -1;
392 Float_t amp = -1;
356fd0a9 393 Float_t timeR = -1;
f0377b23 394 Int_t iSupMod = -1;
395 Int_t nTower = -1;
396 Int_t nIphi = -1;
397 Int_t nIeta = -1;
398 Int_t iphi = -1;
399 Int_t ieta = -1;
356fd0a9 400
401 //List of TRU matrices initialized to 0.
402 for(Int_t k = 0; k < fNTRU*fNumberOfSuperModules; k++){
403 TMatrixD * amptrus = new TMatrixD(nCellsPhi,nCellsEta) ;
404 TMatrixD * timeRtrus = new TMatrixD(nCellsPhi,nCellsEta) ;
405 for(Int_t i = 0; i < nCellsPhi; i++){
406 for(Int_t j = 0; j < nCellsEta; j++){
407 (*amptrus)(i,j) = 0.0;
408 (*timeRtrus)(i,j) = 0.0;
409 }
410 }
411 new((*ampmatrix)[k]) TMatrixD(*amptrus) ;
412 new((*timeRmatrix)[k]) TMatrixD(*timeRtrus) ;
413 }
414
415 AliEMCALDigit * dig ;
f0377b23 416
417 //Digits loop to fill TRU matrices with amplitudes.
f0377b23 418 for(Int_t idig = 0 ; idig < digits->GetEntriesFast() ; idig++){
419
420 dig = dynamic_cast<AliEMCALDigit *>(digits->At(idig)) ;
356fd0a9 421 amp = dig->GetAmp() ; // Energy of the digit (arbitrary units)
422 id = dig->GetId() ; // Id label of the cell
423 timeR = dig->GetTimeR() ; // Earliest time of the digit
f0377b23 424
425 //Get eta and phi cell position in supermodule
426 Bool_t bCell = GetCellIndex(id, iSupMod, nTower, nIphi, nIeta) ;
427 if(!bCell)
428 Error("FillTRU","Wrong cell id number") ;
429
430 GetCellPhiEtaIndexInSModule(iSupMod,nTower,nIphi, nIeta,iphi,ieta);
431
432 //Check to which TRU in the supermodule belongs the cell.
433 //Supermodules are divided in a TRU matrix of dimension
434 //(fNTRUPhi,fNTRUEta).
435 //Each TRU is a cell matrix of dimension (nCellsPhi,nCellsEta)
436
437 //First calculate the row and column in the supermodule
438 //of the TRU to which the cell belongs.
f0377b23 439 Int_t col = (ieta-1)/nCellsEta+1;
440 Int_t row = (iphi-1)/nCellsPhi+1;
356fd0a9 441 if(iSupMod > 10)
442 row = (iphi-1)/nCellsPhi2+1;
443 //Calculate label number of the TRU
444 Int_t itru = (row-1) + (col-1)*fNTRUPhi + (iSupMod-1)*fNTRU ;
f0377b23 445
446 //Fill TRU matrix with cell values
356fd0a9 447 TMatrixD * amptrus = dynamic_cast<TMatrixD *>(ampmatrix->At(itru)) ;
448 TMatrixD * timeRtrus = dynamic_cast<TMatrixD *>(timeRmatrix->At(itru)) ;
f0377b23 449
356fd0a9 450 //Calculate row and column of the cell inside the TRU with number itru
451 Int_t irow = (iphi-1) - (row-1) * nCellsPhi;
452 if(iSupMod > 10)
453 irow = (iphi-1) - (row-1) * nCellsPhi2;
f0377b23 454 Int_t icol = (ieta-1) - (col-1) * nCellsEta;
455
356fd0a9 456 (*amptrus)(irow,icol) = amp ;
457 (*timeRtrus)(irow,icol) = timeR ;
f0377b23 458
459 }
f0377b23 460}
461
356fd0a9 462//______________________________________________________________________
463void AliEMCALGeometry::GetCellPhiEtaIndexInSModuleFromTRUIndex(const Int_t itru, const Int_t iphitru, const Int_t ietatru, Int_t &iphiSM, Int_t &ietaSM) const
464{
465
466 // This method transforms the (eta,phi) index of a cells in a
467 // TRU matrix into Super Module (eta,phi) index.
468
469 // Calculate in which row and column in which the TRU are
470 // ordered in the SM
471
472 Int_t col = itru/ fNTRUPhi + 1;
473 Int_t row = itru - (col-1)*fNTRUPhi + 1;
474
475 //Calculate the (eta,phi) index in SM
476 Int_t nCellsPhi = fNPhi*2/fNTRUPhi;
477 Int_t nCellsEta = fNZ*2/fNTRUEta;
478
479 iphiSM = nCellsPhi*(row-1) + iphitru + 1 ;
480 ietaSM = nCellsEta*(col-1) + ietatru + 1 ;
481}
f0377b23 482
b13bbe81 483//______________________________________________________________________
484AliEMCALGeometry * AliEMCALGeometry::GetInstance(){
05a92d59 485 // Returns the pointer of the unique instance
486
e52475ed 487 AliEMCALGeometry * rv = static_cast<AliEMCALGeometry *>( fgGeom );
488 return rv;
2012850d 489}
173558f2 490
b13bbe81 491//______________________________________________________________________
492AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
493 const Text_t* title){
494 // Returns the pointer of the unique instance
495
496 AliEMCALGeometry * rv = 0;
497 if ( fgGeom == 0 ) {
498 if ( strcmp(name,"") == 0 ) rv = 0;
e5a93224 499 else {
b13bbe81 500 fgGeom = new AliEMCALGeometry(name, title);
501 if ( fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
502 else {
503 rv = 0;
504 delete fgGeom;
505 fgGeom = 0;
506 } // end if fgInit
507 } // end if strcmp(name,"")
508 }else{
e5a93224 509 if ( strcmp(fgGeom->GetName(), name) != 0) {
510 printf("\ncurrent geometry is %s : ", fgGeom->GetName());
511 printf(" you cannot call %s ", name);
b13bbe81 512 }else{
9859bfc0 513 rv = (AliEMCALGeometry *) fgGeom;
e52475ed 514 } // end
b13bbe81 515 } // end if fgGeom
516 return rv;
2012850d 517}
173558f2 518
ab37d09c 519Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const {
dc7da436 520 // Checks whether point is inside the EMCal volume, used in AliEMCALv*.cxx
ab37d09c 521 //
522 // Code uses cylindrical approximation made of inner radius (for speed)
523 //
524 // Points behind EMCAl, i.e. R > outer radius, but eta, phi in acceptance
525 // are considered to inside
526
527 Double_t r=sqrt(x*x+y*y);
528
529 if ( r > fEnvelop[0] ) {
530 Double_t theta;
531 theta = TMath::ATan2(r,z);
532 Double_t eta;
533 if(theta == 0)
534 eta = 9999;
535 else
536 eta = -TMath::Log(TMath::Tan(theta/2.));
537 if (eta < fArm1EtaMin || eta > fArm1EtaMax)
538 return 0;
539
540 Double_t phi = TMath::ATan2(y,x) * 180./TMath::Pi();
541 if (phi > fArm1PhiMin && phi < fArm1PhiMax)
542 return 1;
543 }
544 return 0;
545}
c63c3c5d 546// ==
1963b290 547
548//
549// == Shish-kebab cases ==
550//
e52475ed 551Int_t AliEMCALGeometry::GetAbsCellId(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta) const
dc7da436 552{
553 // 27-aug-04;
d87bd045 554 // corr. 21-sep-04;
555 // 13-oct-05; 110 degree case
dc7da436 556 // May 31, 2006; ALICE numbering scheme:
557 // 0 <= nSupMod < fNumberOfSuperModules
558 // 0 <= nTower < fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1)
559 // 0 <= nIphi < fNPHIdiv
560 // 0 <= nIeta < fNETAdiv
561 // 0 <= absid < fNCells
562 static Int_t id=0; // have to change from 0 to fNCells-1
563 if(fKey110DEG == 1 && nSupMod >= 10) { // 110 degree case; last two supermodules
564 id = fNCellsInSupMod*10 + (fNCellsInSupMod/2)*(nSupMod-10);
d87bd045 565 } else {
dc7da436 566 id = fNCellsInSupMod*nSupMod;
d87bd045 567 }
dc7da436 568 id += fNCellsInTower *nTower;
569 id += fNPHIdiv *nIphi;
1963b290 570 id += nIeta;
dc7da436 571 if(id<0 || id >= fNCells) {
500aeccc 572// printf(" wrong numerations !!\n");
573// printf(" id %6i(will be force to -1)\n", id);
574// printf(" fNCells %6i\n", fNCells);
575// printf(" nSupMod %6i\n", nSupMod);
576// printf(" nTower %6i\n", nTower);
577// printf(" nIphi %6i\n", nIphi);
578// printf(" nIeta %6i\n", nIeta);
dc7da436 579 id = -TMath::Abs(id); // if negative something wrong
1963b290 580 }
581 return id;
582}
583
dc7da436 584Bool_t AliEMCALGeometry::CheckAbsCellId(Int_t absId) const
fc575e27 585{
dc7da436 586 // May 31, 2006; only trd1 now
587 if(absId<0 || absId >= fNCells) return kFALSE;
588 else return kTRUE;
1963b290 589}
590
e52475ed 591Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nTower,Int_t &nIphi,Int_t &nIeta) const
fc575e27 592{
dc7da436 593 // 21-sep-04; 19-oct-05;
594 // May 31, 2006; ALICE numbering scheme:
d87bd045 595 static Int_t tmp=0, sm10=0;
dc7da436 596 if(!CheckAbsCellId(absId)) return kFALSE;
597
d87bd045 598 sm10 = fNCellsInSupMod*10;
dc7da436 599 if(fKey110DEG == 1 && absId >= sm10) { // 110 degree case; last two supermodules
600 nSupMod = (absId-sm10) / (fNCellsInSupMod/2) + 10;
601 tmp = (absId-sm10) % (fNCellsInSupMod/2);
d87bd045 602 } else {
dc7da436 603 nSupMod = absId / fNCellsInSupMod;
604 tmp = absId % fNCellsInSupMod;
d87bd045 605 }
1963b290 606
dc7da436 607 nTower = tmp / fNCellsInTower;
1963b290 608 tmp = tmp % fNCellsInTower;
dc7da436 609 nIphi = tmp / fNPHIdiv;
610 nIeta = tmp % fNPHIdiv;
1963b290 611
612 return kTRUE;
613}
614
dc7da436 615void AliEMCALGeometry::GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, int &iphim, int &ietam) const
fc575e27 616{
dc7da436 617 // added nSupMod; have to check - 19-oct-05 !
618 // Alice numbering scheme - Jun 01,2006
d87bd045 619 static Int_t nphi;
620
dc7da436 621 if(fKey110DEG == 1 && nSupMod>=10) nphi = fNPhi/2;
d87bd045 622 else nphi = fNPhi;
623
dc7da436 624 ietam = nTower/nphi; // have to change from 0 to fNZ-1
625 iphim = nTower%nphi; // have to change from 0 to fNPhi-1
d87bd045 626}
627
eb0b1051 628void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta,
e52475ed 629int &iphi, int &ieta) const
fc575e27 630{
631 // added nSupMod; Nov 25, 05
dc7da436 632 // Alice numbering scheme - Jun 01,2006
633 static Int_t iphim, ietam;
634
635 GetModulePhiEtaIndexInSModule(nSupMod,nTower, iphim, ietam);
636 // have to change from 0 to (fNZ*fNETAdiv-1)
637 ieta = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM)
638 // iphi - have to change from 0 to (fNPhi*fNPHIdiv-1)
639 iphi = iphim*fNPHIdiv + nIphi; // y(module) = y(SM)
1963b290 640}
e52475ed 641
642Int_t AliEMCALGeometry::GetSuperModuleNumber(Int_t absId) const
643{
fc575e27 644 //return the number of the
645 //supermodule given the absolute
646 //ALICE numbering
647
e52475ed 648 static Int_t nSupMod, nTower, nIphi, nIeta;
649 GetCellIndex(absId, nSupMod, nTower, nIphi, nIeta);
650 return nSupMod;
651}
652
653// Methods for AliEMCALRecPoint - Feb 19, 2006
654Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr)
655{
fc575e27 656 //Look to see what the relative
657 //position inside a given cell is
658 //for a recpoint.
659
e52475ed 660 static Int_t nSupMod, nTower, nIphi, nIeta, iphi, ieta;
661 if(!CheckAbsCellId(absId)) return kFALSE;
662
663 GetCellIndex(absId, nSupMod, nTower, nIphi, nIeta);
664 GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi, ieta);
665
dc7da436 666 xr = fXCentersOfCells.At(ieta);
667 zr = fEtaCentersOfCells.At(ieta);
e52475ed 668
dc7da436 669 yr = fPhiCentersOfCells.At(iphi);
e52475ed 670
671 // cout<<" absId "<<absId<<" iphi "<<iphi<<"ieta"<<ieta;
672 // cout<< " xr " << xr << " yr " << yr << " zr " << zr <<endl;
673 return kTRUE;
674}
675
676void AliEMCALGeometry::CreateListOfTrd1Modules()
677{
fc575e27 678 //Generate the list of Trd1 modules
679 //which will make up the EMCAL
680 //geometry
681
e5a93224 682 AliDebug(2,Form(" AliEMCALGeometry::CreateListOfTrd1Modules() started "));
683
e52475ed 684 AliEMCALShishKebabTrd1Module *mod=0, *mTmp=0; // current module
685 if(fShishKebabTrd1Modules == 0) {
686 fShishKebabTrd1Modules = new TList;
687 for(int iz=0; iz< GetNZ(); iz++) {
688 if(iz==0) {
689 mod = new AliEMCALShishKebabTrd1Module(TMath::Pi()/2.,this);
690 } else {
691 mTmp = new AliEMCALShishKebabTrd1Module(*mod);
692 mod = mTmp;
693 }
694 fShishKebabTrd1Modules->Add(mod);
695 }
696 } else {
e5a93224 697 AliDebug(2,Form(" Already exits : "));
e52475ed 698 }
e5a93224 699 AliDebug(2,Form(" fShishKebabTrd1Modules has %i modules \n",
700 fShishKebabTrd1Modules->GetSize()));
e52475ed 701 // Feb 20,2006;
dc7da436 702 // Jun 01, 2006 - ALICE numbering scheme
e52475ed 703 // define grid for cells in eta(z) and x directions in local coordinates system of SM
dc7da436 704 // fEtaCentersOfCells = new TArrayD(fNZ *fNETAdiv);
705 // fXCentersOfCells = new TArrayD(fNZ *fNETAdiv);
706 fEtaCentersOfCells.Set(fNZ *fNETAdiv);
707 fXCentersOfCells.Set(fNZ *fNETAdiv);
708 AliDebug(2,Form(" Cells grid in eta directions : size %i\n", fEtaCentersOfCells.GetSize()));
e52475ed 709 Int_t iphi=0, ieta=0, nTower=0;
710 Double_t xr, zr;
711 for(Int_t it=0; it<fNZ; it++) { // array index
712 AliEMCALShishKebabTrd1Module *trd1 = GetShishKebabModule(it);
dc7da436 713 nTower = fNPhi*it;
e52475ed 714 for(Int_t ic=0; ic<fNETAdiv; ic++) { // array index
dc7da436 715 trd1->GetCenterOfCellInLocalCoordinateofSM(ic, xr, zr);
716 GetCellPhiEtaIndexInSModule(0, nTower, 0, ic, iphi, ieta); // don't depend from phi - ieta in action
717 fXCentersOfCells.AddAt(float(xr) - fParSM[0],ieta);
718 fEtaCentersOfCells.AddAt(float(zr) - fParSM[2],ieta);
e52475ed 719 }
720 }
dc7da436 721 for(Int_t i=0; i<fEtaCentersOfCells.GetSize(); i++) {
e5a93224 722 AliDebug(2,Form(" ind %2.2i : z %8.3f : x %8.3f", i+1,
dc7da436 723 fEtaCentersOfCells.At(i),fXCentersOfCells.At(i)));
e52475ed 724 }
e5a93224 725
e52475ed 726 // define grid for cells in phi(y) direction in local coordinates system of SM
dc7da436 727 // fPhiCentersOfCells = new TArrayD(fNPhi*fNPHIdiv);
728 fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
729 AliDebug(2,Form(" Cells grid in phi directions : size %i\n", fPhiCentersOfCells.GetSize()));
e52475ed 730 Int_t ind=0;
731 for(Int_t it=0; it<fNPhi; it++) { // array index
732 Float_t ytLeftCenterModule = -fParSM[1] + fPhiModuleSize*(2*it+1)/2; // module
733 for(Int_t ic=0; ic<fNPHIdiv; ic++) { // array index
734 Float_t ytLeftCenterCell = ytLeftCenterModule + fPhiTileSize *(2*ic-1)/2.; // tower(cell)
dc7da436 735 fPhiCentersOfCells.AddAt(ytLeftCenterCell,ind);
736 AliDebug(2,Form(" ind %2.2i : y %8.3f ", ind, fPhiCentersOfCells.At(ind)));
e5a93224 737 ind++;
e52475ed 738 }
739 }
e52475ed 740}
741
742void AliEMCALGeometry::GetTransformationForSM()
743{
fc575e27 744 //Uses the geometry manager to
745 //load the transformation matrix
746 //for the supermodules
747
e52475ed 748 static Bool_t transInit=kFALSE;
749 if(transInit) return;
750
751 int i=0;
752 if(gGeoManager == 0) {
753 Info("CreateTransformationForSM() "," Load geometry : TGeoManager::Import()");
754 assert(0);
755 }
756 TGeoNode *tn = gGeoManager->GetTopNode();
fc575e27 757 TGeoNode *node=0, *xen1 = 0;
e52475ed 758 for(i=0; i<tn->GetNdaughters(); i++) {
759 node = tn->GetDaughter(i);
760 TString ns(node->GetName());
761 if(ns.Contains(GetNameOfEMCALEnvelope())) {
fc575e27 762 xen1 = node;
e52475ed 763 break;
764 }
765 }
fc575e27 766 if(!xen1) {
e52475ed 767 Info("CreateTransformationForSM() "," geometry has not EMCAL envelope with name %s",
768 GetNameOfEMCALEnvelope());
769 assert(0);
770 }
fc575e27 771 printf(" i %i : EMCAL Envelope is %s : #SM %i \n", i, xen1->GetName(), xen1->GetNdaughters());
772 for(i=0; i<xen1->GetNdaughters(); i++) {
773 TGeoNodeMatrix *sm = (TGeoNodeMatrix*)xen1->GetDaughter(i);
e52475ed 774 fMatrixOfSM[i] = sm->GetMatrix();
64942713 775 //Compiler doesn't like this syntax...
776 // printf(" %i : matrix %x \n", i, fMatrixOfSM[i]);
e52475ed 777 }
778 transInit = kTRUE;
779}
780
781void AliEMCALGeometry::GetGlobal(const Double_t *loc, Double_t *glob, int nsm) const
782{
fc575e27 783 //Figure out the global numbering
784 //of a given supermodule from the
785 //local numbering
786
e52475ed 787 // if(fMatrixOfSM[0] == 0) GetTransformationForSM();
788 static int ind;
789 ind = nsm-1;
790 if(ind>=0 && ind < GetNumberOfSuperModules()) {
791 fMatrixOfSM[ind]->LocalToMaster(loc, glob);
792 }
793}
794
21cf2b24 795void AliEMCALGeometry::GetGlobal(Int_t /* absId */, TVector3 & /* vglob */) const
e52475ed 796{ // have to be defined
797}
798
799void AliEMCALGeometry::GetGlobal(const TVector3 &vloc, TVector3 &vglob, int nsm) const
800{
fc575e27 801 //Figure out the global numbering
802 //of a given supermodule from the
803 //local numbering given a 3-vector location
804
e52475ed 805 static Double_t tglob[3], tloc[3];
806 vloc.GetXYZ(tloc);
807 GetGlobal(tloc, tglob, nsm);
808 vglob.SetXYZ(tglob[0], tglob[1], tglob[2]);
809}
810
811void AliEMCALGeometry::GetGlobal(const AliRecPoint *rp, TVector3 &vglob) const
812{
664bfd66 813 // Figure out the global numbering
814 // of a given supermodule from the
815 // local numbering for RecPoints
fc575e27 816
e52475ed 817 static TVector3 vloc;
818 static Int_t nSupMod, nTower, nIphi, nIeta;
819
820 AliRecPoint *rpTmp = (AliRecPoint*)rp; // const_cast ??
821 if(!rpTmp) return;
822 AliEMCALRecPoint *rpEmc = (AliEMCALRecPoint*)rpTmp;
823
824 GetCellIndex(rpEmc->GetAbsId(0), nSupMod, nTower, nIphi, nIeta);
825 rpTmp->GetLocalPosition(vloc);
826 GetGlobal(vloc, vglob, nSupMod);
827}
828
664bfd66 829void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Float_t &eta,Float_t &phi) const
830{
831 // Jun 03, 2006 - version for TRD1
832 static TVector3 vglob;
833 GetGlobal(absId, vglob);
834 eta = vglob.Eta();
835 phi = vglob.Phi();
836}
837
fc575e27 838AliEMCALShishKebabTrd1Module* AliEMCALGeometry::GetShishKebabModule(Int_t neta=0)
839{
840 //This method was too long to be
841 //included in the header file - the
842 //rule checker complained about it's
843 //length, so we move it here. It returns the
844 //shishkebabmodule at a given eta index point.
845
846 static AliEMCALShishKebabTrd1Module* trd1=0;
847 if(fShishKebabTrd1Modules && neta>=0 && neta<fShishKebabTrd1Modules->GetSize()) {
848 trd1 = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(neta);
849 } else trd1 = 0;
850 return trd1;
851}