fixed method from id to angle (Bjorn said)
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALGeometry.cxx
1 /**************************************************************************
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  
20 // EMCAL consists of layers of scintillator and lead
21 // Places the the Barrel Geometry of The EMCAL at Midrapidity
22 // between 0 and 120 degrees of Phi and
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               
26 // EMCALArch2x has more modules along both phi and eta
27 // EMCALArchxa has less Layers in the Radial Direction
28 //*-- Author: Sahal Yacoob (LBL / UCT)
29 //     and  : Yves Schutz (SUBATECH)
30 //     and  : Jennifer Klay (LBL)
31
32 // --- ROOT system ---
33
34 // --- Standard library ---
35
36 #include <iostream.h>
37
38 // --- AliRoot header files ---
39 #include <TMath.h>
40 // -- ALICE Headers.
41 #include "AliConst.h"
42 // --- EMCAL headers
43 #include "AliEMCALGeometry.h"
44
45 ClassImp(AliEMCALGeometry);
46
47 AliEMCALGeometry *AliEMCALGeometry::fgGeom = 0;
48 Bool_t            AliEMCALGeometry::fgInit = kFALSE;
49
50 //______________________________________________________________________
51 AliEMCALGeometry::~AliEMCALGeometry(void){
52     // dtor
53 }
54 //______________________________________________________________________
55 void AliEMCALGeometry::Init(void){
56     // Initializes the EMCAL parameters
57
58     fgInit = kFALSE; // Assume failer untill proven otherwise.
59     if(!(  (strcmp( fName, "EMCALArch1a" ) == 0) |
60            (strcmp( fName, "EMCALArch1b" ) == 0) | 
61            (strcmp( fName, "EMCALArch2a" ) == 0) | 
62            (strcmp( fName, "EMCALArch2b" ) == 0) )){
63         cout <<"Instance " << fName << " undefined" << endl;
64     } // end if
65     //
66     if (((strcmp( fName, "EMCALArch1a" ))    == 0) |
67         ((strcmp( fName, "EMCALArch1b" ))    == 0)){
68         fNZ         = 96;
69         fNPhi       = 144;
70     } // end if
71     if (((strcmp( fName, "EMCALArch2a" ))    == 0) |
72         ((strcmp( fName, "EMCALArch2b" ))    == 0)){
73         fNZ         = 112;
74         fNPhi       = 168;
75     } // end if
76     if (((strcmp( fName, "EMCALArch1a" ))    == 0) |
77         ((strcmp( fName, "EMCALArch2a" ))    == 0)){
78         fNLayers    = 21;
79     } // end if
80     if (((strcmp( fName, "EMCALArch1b" ))    == 0) |
81         ((strcmp( fName, "EMCALArch2b" ))    == 0)){
82         fNLayers    = 25;
83     } // end if
84
85     // geometry
86     fAirGap         = 5.0; // cm, air gap between EMCAL mother volume and 
87                            // active material.
88     fAlFrontThick   = 3.18; // cm, Thickness of front Al layer
89     fPbRadThickness = 0.5; // cm, Thickness of theh Pb radiators.
90     fPreShowerSintThick = 0.6; // cm, Thickness of the sintilator for the
91                                // preshower part of the calorimeter
92     fFullShowerSintThick = 0.5; // cm, Thickness of the sintilator for the
93                                 // full shower part of the calorimeter
94     fArm1PhiMin     =   0.0; // degrees, Starting EMCAL Phi position
95     fArm1PhiMax     = 120.0; // degrees, Ending EMCAL Phi position
96     fArm1EtaMin     = -0.7; // pseudorapidity, Starting EMCAL Eta position
97     fArm1EtaMax     = +0.7; // pseudorapidity, Ending EMCAL Eta position
98     fIPDistance     = 454.0; // cm, Radial distance to inner surface of EMCAL
99     fShellThickness = GetAlFrontThickness() + 2.*GetPreSintThick() +
100         (fNLayers-2)*GetFullSintThick()+(fNLayers-1)*GetPbRadThick();
101     //below; cm, Z lenght of the EMCAL.
102     fZLength        = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax);
103     fEnvelop[0]     = fIPDistance; // mother volume inner radius
104     fEnvelop[1]     = fIPDistance + fShellThickness; // mother volume outer r.
105     fEnvelop[2]     = 1.00001*fZLength; // add some padding for mother volume. 
106     fGap2Active     = 1.0;  // cm, Gap between 
107     fgInit = kTRUE; 
108 }
109 //______________________________________________________________________
110 AliEMCALGeometry *  AliEMCALGeometry::GetInstance(){ 
111     // Returns the pointer of the unique instance
112
113     return (AliEMCALGeometry *) fgGeom; 
114 }
115 //______________________________________________________________________
116 AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
117                                                 const Text_t* title){
118     // Returns the pointer of the unique instance
119
120     AliEMCALGeometry * rv = 0; 
121     if ( fgGeom == 0 ) {
122         if ( strcmp(name,"") == 0 ) rv = 0;
123         else {    
124             fgGeom = new AliEMCALGeometry(name, title);
125             if ( fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
126             else {
127                 rv = 0; 
128                 delete fgGeom; 
129                 fgGeom = 0; 
130             } // end if fgInit
131         } // end if strcmp(name,"")
132     }else{
133         if ( strcmp(fgGeom->GetName(), name) != 0 ) {
134             cout << "AliEMCALGeometry <E> : current geometry is " 
135                  << fgGeom->GetName() << endl
136                  << "                      you cannot call     " << name 
137                  << endl; 
138         }else{
139             rv = (AliEMCALGeometry *) fgGeom; 
140         } // end if
141     }  // end if fgGeom
142     return rv; 
143 }
144 //______________________________________________________________________
145 Int_t AliEMCALGeometry::TowerIndex(Int_t ieta,Int_t iphi,Int_t ipre){
146     // Returns the tower index number from the based on the Z and Phi
147     // index numbers. There are 2 times the number of towers to separate
148     // out the full towsers from the pre-towsers.
149     // Inputs:
150     //   Int_t ieta    // index allong z axis [1-fNZ]
151     //   Int_t iphi  // index allong phi axis [1-fNPhi]
152     //   Int_t ipre  // 0 = Full tower, 1 = Pre-shower tower only. [0,1]
153     // Outputs:
154     //   none.
155     // Returned
156     // Int_t the absoulute tower index. [1-2*fNZ*fNPhi]
157     Int_t index;
158
159     if((ieta<=0 || ieta>GetNEta()) || (iphi<=0 || iphi>GetNPhi()) ||
160        (ipre<0 || ipre>1) ){
161         cout << "inputs out of range ieta=" << ieta << " [1-" << GetNEta();
162         cout << "] iphi=" << iphi << " [1-" << GetNPhi() << "] ipre=";
163         cout << ipre << "[0,1]. returning -1" << endl;
164         return -1;
165     } // end if
166     index = iphi + GetNPhi()*(ieta-1) + ipre*(GetNPhi()*GetNEta());
167     return index;
168 }
169 //______________________________________________________________________
170 void AliEMCALGeometry::TowerIndexes(Int_t index,Int_t &ieta,Int_t &iphi,
171                                     Int_t &ipre){
172     // given the tower index number it returns the based on the Z and Phi
173     // index numbers and if it is for the full tower or the pre-tower number.
174     // There are 2 times the number of towers to separate
175     // out the full towsers from the pre-towsers.
176     // Inputs:
177     //   Int_t index // Tower index number [1-2*fNZ*fNPhi]
178     // Outputs:
179     //   Int_t ieta    // index allong z axis [1-fNZ]
180     //   Int_t iphi  // index allong phi axis [1-fNPhi]
181     //   Int_t ipre  // 0 = Full tower, 1 = Pre-shower tower only. [0,1]
182     // Returned
183     //   none.
184     Int_t itowers;
185
186     itowers = GetNEta()*GetNPhi();
187     if(index<1 || index>2*itowers){
188         cout << "index=" << index <<" is out of range [1-";
189         cout << 2*itowers << "], returning -1 for all." << endl;
190         ieta = -1; iphi = -1; ipre = -1;
191         return ;
192     } // end if
193     ipre = 0;
194     if(index>itowers){ // pre shower indexs
195         ipre = 1;
196         index = index - itowers;
197     } // end if
198     ieta = 1+ (Int_t)((index-1)/GetNPhi());
199     iphi = index - GetNPhi()*(ieta-1);
200     return;
201 }
202 //______________________________________________________________________
203 void AliEMCALGeometry::EtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi){
204     // given the tower index number it returns the based on the eta and phi
205     // of the tower.
206     // Inputs:
207     //   Int_t index // Tower index number [1-2*fNZ*fNPhi]
208     // Outputs:
209     //   Float_t eta  // eta of center of tower in pseudorapidity
210     //   Float_t phi  // phi of center of tower in degrees
211     // Returned
212     //   none.
213     Int_t ieta,iphi,ipre;
214     Double_t deta,dphi,phid;
215
216     TowerIndexes(index,ieta,iphi,ipre);
217     deta = (GetArm1EtaMax()-GetArm1EtaMin())/((Float_t)GetNEta());
218     eta  = GetArm1EtaMin() + (((Float_t)ieta)-0.5)*deta;
219     dphi = (GetArm1PhiMax() - GetArm1PhiMin())/((Float_t)GetNPhi());  // in degrees.
220     phid = GetArm1PhiMin() + dphi*((Float_t)iphi -0.5);//iphi range [1-fNphi].
221     phi  = phid;
222 }
223 //______________________________________________________________________
224 Int_t AliEMCALGeometry::TowerIndexFromEtaPhi(Float_t eta,Float_t phi){
225     // returns the tower index number based on the eta and phi of the tower.
226     // Inputs:
227     //   Float_t eta  // eta of center of tower in pseudorapidity
228     //   Float_t phi  // phi of center of tower in degrees
229     // Outputs:
230     //   none.
231     // Returned
232     //   Int_t index // Tower index number [1-fNZ*fNPhi]
233     Int_t ieta,iphi;
234
235     ieta = 1 + (Int_t)(((Float_t)GetNEta())*(eta-GetArm1EtaMin())/
236                   (GetArm1EtaMax() - GetArm1EtaMin()));
237     if(ieta<=0 || ieta>GetNEta()){
238         cout << "TowerIndexFromEtaPhi:";
239         cout << "ieta = "<< ieta << " eta=" << eta << " is outside of EMCAL. etamin=";
240         cout << GetArm1EtaMin() << " to etamax=" << GetArm1EtaMax();
241         cout << " returning -1" << endl;
242         return -1;
243     } // end if
244     iphi = 1 + (Int_t)(((Float_t)GetNPhi())*(phi-GetArm1PhiMin())/
245                   ((Float_t)(GetArm1PhiMax() - GetArm1PhiMin())));
246     if(iphi<=0 || iphi>GetNPhi()){
247         cout << "TowerIndexFromEtaPhi:";
248         cout << "iphi=" << iphi << " phi=" << phi << " is outside of EMCAL.";
249         cout << " Phimin=" << GetArm1PhiMin() << " PhiMax=" << GetArm1PhiMax();
250         cout << " returning -1" << endl;
251         return -1;
252     } // end if
253     return TowerIndex(ieta,iphi,0);
254 }
255 //______________________________________________________________________
256 Int_t AliEMCALGeometry::PreTowerIndexFromEtaPhi(Float_t eta,Float_t phi){
257     // returns the pretower index number based on the eta and phi of the tower.
258     // Inputs:
259     //   Float_t eta  // eta of center of tower in pseudorapidity
260     //   Float_t phi  // phi of center of tower in degrees
261     // Outputs:
262     //   none.
263     // Returned
264     //   Int_t index // PreTower index number [fNZ*fNPhi-2*fNZ*fNPhi]
265
266     return GetNEta()*GetNPhi()+TowerIndexFromEtaPhi(eta,phi);
267 }
268 //______________________________________________________________________
269 Bool_t AliEMCALGeometry::AbsToRelNumbering(Int_t AbsId, Int_t *relid){
270     // Converts the absolute numbering into the following array/
271     //  relid[0] = EMCAL Arm number 1:1 
272     //  relid[1] = 0  Not in Pre Shower layers
273     //           = -1 In Pre Shower
274     //  relid[2] = Row number inside EMCAL
275     //  relid[3] = Column number inside EMCAL
276     // Input:
277     //   Int_t AbsId // Tower index number [1-2*fNZ*fNPhi]
278     // Outputs:
279     //   Int_t *relid // array of 5. Discribed above.
280     Bool_t rv  = kTRUE ;
281     Int_t ieta=0,iphi=0,ipre=0,index=AbsId;
282
283     TowerIndexes(index,ieta,iphi,ipre);
284     relid[0] = 1;
285     relid[1] = 0;
286     if(ipre==1) 
287       relid[1] = -1;
288     relid[2] = ieta;
289     relid[3] = iphi;
290
291     return rv;
292 }
293 //______________________________________________________________________
294 void AliEMCALGeometry::PosInAlice(const Int_t *relid,Float_t &theta,
295                                      Float_t &phi){
296     // Converts the relative numbering into the local EMCAL-module (x, z)
297     // coordinates
298     Int_t ieta   = relid[2]; // offset along x axis
299     Int_t iphi = relid[3]; // offset along z axis
300     Int_t ipre = relid[1]; // indicates -1 preshower, or 0 full tower.
301     Int_t index;
302     Float_t eta;
303
304     if(ipre==-1) ipre = 1;
305     index = TowerIndex(ieta,iphi,ipre);
306     EtaPhiFromIndex(index,eta,phi);
307     theta = 180.*(2.0*TMath::ATan(TMath::Exp(-eta)))/TMath::Pi();
308
309     return;
310 }
311 //______________________________________________________________________
312 /*
313 Boot_t AliEMCALGeometry::AreNeighbours(Int_t index1,Int_t index2){
314     // Returns kTRUE if the two towers are neighbours or not, including
315     // diagonals. Both indexes are required to be either towers or preshower.
316     // Inputs:
317     //   Int_t index1  // index of tower 1
318     //   Int_t index2  // index of tower 2
319     // Outputs:
320     //   none.
321     // Returned
322     //   Boot_t kTRUE if the towers are neighbours otherwise false.
323     Boot_t anb = kFALSE;
324     Int_t ieta1 = 0, ieta2 = 0, iphi1 = 0, iphi2 = 0, ipre1 = 0, ipre2 = 0;
325
326     TowerIndexes(index1,ieta1,iphi1,ipre1);
327     TowerIndexes(index2,ieta2,iphi2,ipre2);
328     if(ipre1!=ipre2) return anb;
329     if((ieta1>=ieta2-1 && ieta1<=ieta2+1) && (iphi1>=iphi2-1 &&iphi1<=iphi2+1))
330                                                                  anb = kTRUE;
331     return anb;
332 }
333  */