removed iostream
[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 #include <stdlib.h> 
36
37 // --- AliRoot header files ---
38 #include <TMath.h>
39
40 // -- ALICE Headers.
41 #include "AliConst.h"
42
43 // --- EMCAL headers
44 #include "AliEMCALGeometry.h"
45
46 ClassImp(AliEMCALGeometry);
47
48 AliEMCALGeometry *AliEMCALGeometry::fgGeom = 0;
49 Bool_t            AliEMCALGeometry::fgInit = kFALSE;
50
51 //______________________________________________________________________
52 AliEMCALGeometry::~AliEMCALGeometry(void){
53     // dtor
54 }
55 //______________________________________________________________________
56 void AliEMCALGeometry::Init(void){
57     // Initializes the EMCAL parameters
58
59     fgInit = kFALSE; // Assume failer untill proven otherwise.
60
61     TString name(GetName()) ; 
62                  
63     if( name != "EMCALArch1a" &&
64         name != "EMCALArch1b" && 
65         name != "EMCALArch2a" && 
66         name != "EMCALArch2b"  ){
67       Fatal("Init", "%s is not a known geometry (choose among EMCALArch1a, EMCALArch1b, EMCALArch2a and EMCALArch2b)",  name.Data()) ;  
68     } // end if
69     //
70     if ( name == "EMCALArch1a"  ||
71          name == "EMCALArch1b" ) {
72         fNZ         = 96;
73         fNPhi       = 144;
74     } // end if
75     if ( name == "EMCALArch2a"  ||
76          name, "EMCALArch2b" ) {
77         fNZ         = 112;
78         fNPhi       = 168;
79     } // end if
80     if ( name == "EMCALArch1a"  ||
81          name == "EMCALArch2a" ) {
82         fNLayers    = 21;
83     } // end if
84     if ( name == "EMCALArch1b"  ||
85          name == "EMCALArch2b" ) {
86         fNLayers    = 25;
87     } // end if
88
89     // geometry
90     fAirGap         = 5.0; // cm, air gap between EMCAL mother volume and 
91                            // active material.
92     fAlFrontThick   = 3.18; // cm, Thickness of front Al layer
93     fPbRadThickness = 0.5; // cm, Thickness of theh Pb radiators.
94     fPreShowerSintThick = 0.6; // cm, Thickness of the sintilator for the
95                                // preshower part of the calorimeter
96     fFullShowerSintThick = 0.5; // cm, Thickness of the sintilator for the
97                                 // full shower part of the calorimeter
98     fArm1PhiMin     =  60.0; // degrees, Starting EMCAL Phi position
99     fArm1PhiMax     = 180.0; // degrees, Ending EMCAL Phi position
100     fArm1EtaMin     = -0.7; // pseudorapidity, Starting EMCAL Eta position
101     fArm1EtaMax     = +0.7; // pseudorapidity, Ending EMCAL Eta position
102     fIPDistance     = 454.0; // cm, Radial distance to inner surface of EMCAL
103     fShellThickness = GetAlFrontThickness() + 2.*GetPreSintThick() +
104         (fNLayers-2)*GetFullSintThick()+(fNLayers-1)*GetPbRadThick();
105     //below; cm, Z lenght of the EMCAL.
106     fZLength        = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax);
107     fEnvelop[0]     = fIPDistance; // mother volume inner radius
108     fEnvelop[1]     = fIPDistance + fShellThickness; // mother volume outer r.
109     fEnvelop[2]     = 1.00001*fZLength; // add some padding for mother volume. 
110     fGap2Active     = 1.0;  // cm, Gap between 
111     fgInit = kTRUE; 
112 }
113
114 //______________________________________________________________________
115 AliEMCALGeometry *  AliEMCALGeometry::GetInstance(){ 
116   // Returns the pointer of the unique instance
117   
118   return static_cast<AliEMCALGeometry *>( fgGeom ) ; 
119 }
120
121 //______________________________________________________________________
122 AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
123                                                 const Text_t* title){
124     // Returns the pointer of the unique instance
125
126     AliEMCALGeometry * rv = 0; 
127     if ( fgGeom == 0 ) {
128         if ( strcmp(name,"") == 0 ) rv = 0;
129         else {    
130             fgGeom = new AliEMCALGeometry(name, title);
131             if ( fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
132             else {
133                 rv = 0; 
134                 delete fgGeom; 
135                 fgGeom = 0; 
136             } // end if fgInit
137         } // end if strcmp(name,"")
138     }else{
139         if ( strcmp(fgGeom->GetName(), name) != 0 ) {
140           TString message("\n") ; 
141           message += "current geometry is " ;  
142           message += fgGeom->GetName() ;
143           message += "\n                      you cannot call     " ; 
144           message += name ;  
145           ::Info("GetGeometry", message.Data() ) ; 
146         }else{
147           rv = (AliEMCALGeometry *) fgGeom; 
148         } // end if
149     }  // end if fgGeom
150     return rv; 
151 }
152
153 //______________________________________________________________________
154 Int_t AliEMCALGeometry::TowerIndex(Int_t ieta,Int_t iphi,Int_t ipre) const {
155     // Returns the tower index number from the based on the Z and Phi
156     // index numbers. There are 2 times the number of towers to separate
157     // out the full towsers from the pre-towsers.
158     // Inputs:
159     //   Int_t ieta    // index allong z axis [1-fNZ]
160     //   Int_t iphi  // index allong phi axis [1-fNPhi]
161     //   Int_t ipre  // 0 = Full tower, 1 = Pre-shower tower only. [0,1]
162     // Outputs:
163     //   none.
164     // Returned
165     // Int_t the absoulute tower index. [1-2*fNZ*fNPhi]
166     Int_t index;
167
168     if((ieta<=0 || ieta>GetNEta()) || (iphi<=0 || iphi>GetNPhi()) ||
169        (ipre<0 || ipre>1) ){
170       TString message ("\n") ; 
171       message += "inputs out of range ieta= " ; 
172       message += ieta ; 
173       message += " [1-" ; 
174       message += GetNEta() ;
175       message += "] iphi= " ; 
176       message += iphi ; 
177       message += " [1-" ; 
178       message += GetNPhi() ; 
179       message += "] ipre= " ;
180       message += ipre ; 
181       message += "[0,1]. returning -1" ; 
182       Warning("TowerIndex", message.Data() ) ; 
183       return -1;
184     } // end if
185     index = iphi + GetNPhi()*(ieta-1) + ipre*(GetNPhi()*GetNEta());
186     return index;
187 }
188
189 //______________________________________________________________________
190 void AliEMCALGeometry::TowerIndexes(Int_t index,Int_t &ieta,Int_t &iphi,
191                                     Int_t &ipre) const {
192     // given the tower index number it returns the based on the Z and Phi
193     // index numbers and if it is for the full tower or the pre-tower number.
194     // There are 2 times the number of towers to separate
195     // out the full towsers from the pre-towsers.
196     // Inputs:
197     //   Int_t index // Tower index number [1-2*fNZ*fNPhi]
198     // Outputs:
199     //   Int_t ieta    // index allong z axis [1-fNZ]
200     //   Int_t iphi  // index allong phi axis [1-fNPhi]
201     //   Int_t ipre  // 0 = Full tower, 1 = Pre-shower tower only. [0,1]
202     // Returned
203     //   none.
204     Int_t itowers;
205
206     itowers = GetNEta()*GetNPhi();
207     if(index<1 || index>2*itowers){
208       TString message("\n") ; 
209       message += "index= " ; 
210       message += index ; 
211       message += " is out of range [1-" ;
212       message += 2*itowers ; 
213       message += "], returning -1 for all." ;
214       Warning("TowerIndex", message.Data() ) ; 
215       ieta = -1; iphi = -1; ipre = -1;
216       return ;
217     } // end if
218     ipre = 0;
219     if(index>itowers){ // pre shower indexs
220         ipre = 1;
221         index = index - itowers;
222     } // end if
223     ieta = 1+ (Int_t)((index-1)/GetNPhi());
224     iphi = index - GetNPhi()*(ieta-1);
225     return;
226 }
227
228 //______________________________________________________________________
229 void AliEMCALGeometry::EtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi) const {
230     // given the tower index number it returns the based on the eta and phi
231     // of the tower.
232     // Inputs:
233     //   Int_t index // Tower index number [1-2*fNZ*fNPhi]
234     // Outputs:
235     //   Float_t eta  // eta of center of tower in pseudorapidity
236     //   Float_t phi  // phi of center of tower in degrees
237     // Returned
238     //   none.
239     Int_t ieta,iphi,ipre;
240     Double_t deta,dphi,phid;
241
242     TowerIndexes(index,ieta,iphi,ipre);
243     deta = (GetArm1EtaMax()-GetArm1EtaMin())/((Float_t)GetNEta());
244     eta  = GetArm1EtaMin() + (((Float_t)ieta)-0.5)*deta;
245     dphi = (GetArm1PhiMax() - GetArm1PhiMin())/((Float_t)GetNPhi());  // in degrees.
246     phid = GetArm1PhiMin() + dphi*((Float_t)iphi -0.5);//iphi range [1-fNphi].
247     phi  = phid;
248 }
249
250 //______________________________________________________________________
251 Int_t AliEMCALGeometry::TowerIndexFromEtaPhi(Float_t eta,Float_t phi) const {
252     // returns the tower index number based on the eta and phi of the tower.
253     // Inputs:
254     //   Float_t eta  // eta of center of tower in pseudorapidity
255     //   Float_t phi  // phi of center of tower in degrees
256     // Outputs:
257     //   none.
258     // Returned
259     //   Int_t index // Tower index number [1-fNZ*fNPhi]
260     Int_t ieta,iphi;
261
262     ieta = 1 + (Int_t)(((Float_t)GetNEta())*(eta-GetArm1EtaMin())/
263                   (GetArm1EtaMax() - GetArm1EtaMin()));
264     if(ieta<=0 || ieta>GetNEta()){
265       TString message("\n") ; 
266       message += "ieta = " ; 
267       message += ieta ; 
268       message += " eta=" ; 
269       message += eta ; 
270       message += " is outside of EMCAL. etamin=" ;
271       message += GetArm1EtaMin() ;
272       message += " to etamax=" ; 
273       message += GetArm1EtaMax();
274       message += " returning -1";
275       Warning("TowerIndexFromEtaPhi", message.Data() ) ; 
276       return -1;
277     } // end if
278     iphi = 1 + (Int_t)(((Float_t)GetNPhi())*(phi-GetArm1PhiMin())/
279                   ((Float_t)(GetArm1PhiMax() - GetArm1PhiMin())));
280     if(iphi<=0 || iphi>GetNPhi()){
281       TString message("\n") ; 
282       message += "iphi=" ; 
283       message += iphi ;  
284       message += "phi= " ; 
285       message += phi ; 
286       message += " is outside of EMCAL." ;
287       message += " Phimin=" ; 
288       message += GetArm1PhiMin() ; 
289       message += " PhiMax=" ; 
290       message += GetArm1PhiMax() ;
291       message += " returning -1" ;
292       Warning("TowerIndexFromEtaPhi", message.Data() ) ; 
293       return -1;
294     } // end if
295     return TowerIndex(ieta,iphi,0);
296 }
297
298 //______________________________________________________________________
299 Int_t AliEMCALGeometry::PreTowerIndexFromEtaPhi(Float_t eta,Float_t phi) const {
300     // returns the pretower index number based on the eta and phi of the tower.
301     // Inputs:
302     //   Float_t eta  // eta of center of tower in pseudorapidity
303     //   Float_t phi  // phi of center of tower in degrees
304     // Outputs:
305     //   none.
306     // Returned
307     //   Int_t index // PreTower index number [fNZ*fNPhi-2*fNZ*fNPhi]
308
309     return GetNEta()*GetNPhi()+TowerIndexFromEtaPhi(eta,phi);
310 }
311
312 //______________________________________________________________________
313 Bool_t AliEMCALGeometry::AbsToRelNumbering(Int_t AbsId, Int_t *relid) const {
314     // Converts the absolute numbering into the following array/
315     //  relid[0] = EMCAL Arm number 1:1 
316     //  relid[1] = 0  Not in Pre Shower layers
317     //           = -1 In Pre Shower
318     //  relid[2] = Row number inside EMCAL
319     //  relid[3] = Column number inside EMCAL
320     // Input:
321     //   Int_t AbsId // Tower index number [1-2*fNZ*fNPhi]
322     // Outputs:
323     //   Int_t *relid // array of 5. Discribed above.
324     Bool_t rv  = kTRUE ;
325     Int_t ieta=0,iphi=0,ipre=0,index=AbsId;
326
327     TowerIndexes(index,ieta,iphi,ipre);
328     relid[0] = 1;
329     relid[1] = 0;
330     if(ipre==1) 
331       relid[1] = -1;
332     relid[2] = ieta;
333     relid[3] = iphi;
334
335     return rv;
336 }
337
338 //______________________________________________________________________
339 void AliEMCALGeometry::PosInAlice(const Int_t *relid,Float_t &theta,
340                                      Float_t &phi) const {
341     // Converts the relative numbering into the local EMCAL-module (x, z)
342     // coordinates
343     Int_t ieta   = relid[2]; // offset along x axis
344     Int_t iphi = relid[3]; // offset along z axis
345     Int_t ipre = relid[1]; // indicates -1 preshower, or 0 full tower.
346     Int_t index;
347     Float_t eta;
348
349     if(ipre==-1) ipre = 1;
350     index = TowerIndex(ieta,iphi,ipre);
351     EtaPhiFromIndex(index,eta,phi);
352     theta = 180.*(2.0*TMath::ATan(TMath::Exp(-eta)))/TMath::Pi();
353
354     return;
355 }
356
357 //______________________________________________________________________
358 void AliEMCALGeometry::XYZFromIndex(const Int_t *relid,Float_t &x,Float_t &y, Float_t &z) const {
359     // given the tower relative number it returns the X, Y and Z
360     // of the tower.
361     
362     // Outputs:
363     //   Float_t x  // x of center of tower in cm
364     //   Float_t y  // y of center of tower in cm
365     //   Float_t z  // z of centre of tower in cm
366     // Returned
367     //   none.
368     
369     Float_t eta,theta, phi,cyl_radius,kDeg2Rad;
370     
371     Int_t ieta   = relid[2]; // offset along x axis
372     Int_t iphi = relid[3]; // offset along z axis
373     Int_t ipre = relid[1]; // indicates -1 preshower, or 0 full tower.
374     Int_t index;
375     
376
377     if(ipre==-1) ipre = 1;
378     index = TowerIndex(ieta,iphi,ipre);
379     EtaPhiFromIndex(index,eta,phi);
380     theta = 180.*(2.0*TMath::ATan(TMath::Exp(-eta)))/TMath::Pi();
381     
382     kDeg2Rad = TMath::Pi() / static_cast<Double_t>(180) ; 
383     cyl_radius = GetIPDistance()+ GetAirGap() ;
384     x =  cyl_radius * TMath::Cos(phi * kDeg2Rad ) ;
385     y =  cyl_radius * TMath::Cos(phi * kDeg2Rad ) ; 
386     z =  cyl_radius / TMath::Tan(theta * kDeg2Rad ) ; 
387  
388  return;
389
390
391 //______________________________________________________________________
392 /*
393 Boot_t AliEMCALGeometry::AreNeighbours(Int_t index1,Int_t index2) const {
394     // Returns kTRUE if the two towers are neighbours or not, including
395     // diagonals. Both indexes are required to be either towers or preshower.
396     // Inputs:
397     //   Int_t index1  // index of tower 1
398     //   Int_t index2  // index of tower 2
399     // Outputs:
400     //   none.
401     // Returned
402     //   Boot_t kTRUE if the towers are neighbours otherwise false.
403     Boot_t anb = kFALSE;
404     Int_t ieta1 = 0, ieta2 = 0, iphi1 = 0, iphi2 = 0, ipre1 = 0, ipre2 = 0;
405
406     TowerIndexes(index1,ieta1,iphi1,ipre1);
407     TowerIndexes(index2,ieta2,iphi2,ipre2);
408     if(ipre1!=ipre2) return anb;
409     if((ieta1>=ieta2-1 && ieta1<=ieta2+1) && (iphi1>=iphi2-1 &&iphi1<=iphi2+1))
410                                                                  anb = kTRUE;
411     return anb;
412 }
413  */