Adaption to new fluka common blocks (E. Futo)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSGeometry.cxx
CommitLineData
d15a28e7 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
b2a60966 16/* $Id$ */
17
d15a28e7 18//_________________________________________________________________________
b2a60966 19// Geometry class for PHOS : singleton
a3dfe79c 20// PHOS consists of the electromagnetic calorimeter (EMCA)
21// and a charged particle veto either in the Subatech's version (PPSD)
22// or in the IHEP's one (CPV).
23// The EMCA/PPSD/CPV modules are parametrized so that any configuration
24// can be easily implemented
25// The title is used to identify the version of CPV used.
b2a60966 26//
710f859a 27//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (RRC "KI" & SUBATECH)
d15a28e7 28
29// --- ROOT system ---
30
31#include "TVector3.h"
32#include "TRotation.h"
fa7cce36 33#include "TFolder.h"
34#include "TROOT.h"
d15a28e7 35
36// --- Standard library ---
37
4410223b 38#include <stdlib.h>
d15a28e7 39
40// --- AliRoot header files ---
41
42#include "AliPHOSGeometry.h"
468794ea 43#include "AliPHOSEMCAGeometry.h"
710f859a 44#include "AliPHOSRecPoint.h"
d15a28e7 45#include "AliConst.h"
46
9ec91567 47ClassImp(AliPHOSGeometry) ;
d15a28e7 48
a4e98857 49// these initialisations are needed for a singleton
9ec91567 50AliPHOSGeometry * AliPHOSGeometry::fgGeom = 0 ;
282c5906 51Bool_t AliPHOSGeometry::fgInit = kFALSE ;
9ec91567 52
d15a28e7 53//____________________________________________________________________________
0bc3b8ed 54AliPHOSGeometry::AliPHOSGeometry(void)
55{
56 // default ctor
57 // must be kept public for root persistency purposes,
58 // but should never be called by the outside world
59 fPHOSAngle = 0 ;
60 fGeometryEMCA = 0;
61 fGeometrySUPP = 0;
62 fGeometryCPV = 0;
63 fgGeom = 0;
64 fRotMatrixArray = 0;
65}
66//____________________________________________________________________________
d15a28e7 67AliPHOSGeometry::~AliPHOSGeometry(void)
68{
b2a60966 69 // dtor
70
52a36ffd 71 if (fRotMatrixArray) fRotMatrixArray->Delete() ;
72 if (fRotMatrixArray) delete fRotMatrixArray ;
fa0bc588 73 if (fPHOSAngle ) delete[] fPHOSAngle ;
52a36ffd 74}
52a36ffd 75//____________________________________________________________________________
76
77void AliPHOSGeometry::Init(void)
78{
a4e98857 79 // Initializes the PHOS parameters :
80 // IHEP is the Protvino CPV (cathode pad chambers)
81 // GPS2 is the Subatech Pre-Shower (two micromegas sandwiching a passive lead converter)
809cd394 82 // MIXT 4 PHOS modules withe the IHEP CPV and one PHOS module with the Subatech Pre-Shower
710f859a 83
809cd394 84 TString test(GetName()) ;
85 if (test != "IHEP" && test != "GPS2" && test != "MIXT") {
f1611b7c 86 Fatal("Init", "%s is not a known geometry (choose among IHEP, GPS2 and MIXT)", test.Data() ) ;
809cd394 87 }
88
710f859a 89 fgInit = kTRUE ;
90
91 fNModules = 5;
92 fAngle = 20;
93
94 fGeometryEMCA = new AliPHOSEMCAGeometry();
95
96 fGeometryCPV = new AliPHOSCPVGeometry ();
97
98 fGeometrySUPP = new AliPHOSSupportGeometry();
99
100 fPHOSAngle = new Float_t[fNModules] ;
101
102 Float_t * emcParams = fGeometryEMCA->GetEMCParams() ;
103
581e32d4 104 fPHOSParams[0] = TMath::Max((Double_t)fGeometryCPV->GetCPVBoxSize(0)/2.,
105 (Double_t)(emcParams[0]*(fGeometryCPV->GetCPVBoxSize(1)+emcParams[3]) -
710f859a 106 emcParams[1]* fGeometryCPV->GetCPVBoxSize(1))/emcParams[3] ) ;
107 fPHOSParams[1] = emcParams[1] ;
581e32d4 108 fPHOSParams[2] = TMath::Max((Double_t)emcParams[2], (Double_t)fGeometryCPV->GetCPVBoxSize(2)/2.);
710f859a 109 fPHOSParams[3] = emcParams[3] + fGeometryCPV->GetCPVBoxSize(1)/2. ;
110
111 fIPtoUpperCPVsurface = fGeometryEMCA->GetIPtoOuterCoverDistance() - fGeometryCPV->GetCPVBoxSize(1) ;
112
113 Int_t index ;
114 for ( index = 0; index < fNModules; index++ )
52a36ffd 115 fPHOSAngle[index] = 0.0 ; // Module position angles are set in CreateGeometry()
710f859a 116
117 this->SetPHOSAngles() ;
118 fRotMatrixArray = new TObjArray(fNModules) ;
119
52a36ffd 120}
121
5ccb0008 122//____________________________________________________________________________
52a36ffd 123AliPHOSGeometry * AliPHOSGeometry::GetInstance()
124{
a4e98857 125 // Returns the pointer of the unique instance; singleton specific
126
809cd394 127 return static_cast<AliPHOSGeometry *>( fgGeom ) ;
52a36ffd 128}
129
130//____________________________________________________________________________
131AliPHOSGeometry * AliPHOSGeometry::GetInstance(const Text_t* name, const Text_t* title)
132{
133 // Returns the pointer of the unique instance
a4e98857 134 // Creates it with the specified options (name, title) if it does not exist yet
135
52a36ffd 136 AliPHOSGeometry * rv = 0 ;
137 if ( fgGeom == 0 ) {
138 if ( strcmp(name,"") == 0 )
139 rv = 0 ;
140 else {
141 fgGeom = new AliPHOSGeometry(name, title) ;
142 if ( fgInit )
143 rv = (AliPHOSGeometry * ) fgGeom ;
144 else {
145 rv = 0 ;
146 delete fgGeom ;
147 fgGeom = 0 ;
148 }
149 }
150 }
151 else {
21cd0c07 152 if ( strcmp(fgGeom->GetName(), name) != 0 )
153 ::Error("GetInstance", "Current geometry is %s. You cannot call %s", fgGeom->GetName(), name) ;
52a36ffd 154 else
155 rv = (AliPHOSGeometry *) fgGeom ;
156 }
157 return rv ;
158}
4697edca 159
52a36ffd 160//____________________________________________________________________________
161void AliPHOSGeometry::SetPHOSAngles()
162{
a4e98857 163 // Calculates the position of the PHOS modules in ALICE global coordinate system
52a36ffd 164
165 Double_t const kRADDEG = 180.0 / kPI ;
710f859a 166 Float_t pphi = 2 * TMath::ATan( GetOuterBoxSize(0) / ( 2.0 * GetIPtoUpperCPVsurface() ) ) ;
52a36ffd 167 pphi *= kRADDEG ;
710f859a 168 if (pphi > fAngle){
21cd0c07 169 Error("SetPHOSAngles", "PHOS modules overlap!\n pphi = %f fAngle = %f", pphi, fAngle);
710f859a 170
171 }
ed4205d8 172 pphi = fAngle;
52a36ffd 173
174 for( Int_t i = 1; i <= fNModules ; i++ ) {
ed4205d8 175 Float_t angle = pphi * ( i - fNModules / 2.0 - 0.5 ) ;
52a36ffd 176 fPHOSAngle[i-1] = - angle ;
177 }
d15a28e7 178}
179
180//____________________________________________________________________________
7b7c1533 181Bool_t AliPHOSGeometry::AbsToRelNumbering(const Int_t AbsId, Int_t * relid) const
d15a28e7 182{
b2a60966 183 // Converts the absolute numbering into the following array/
184 // relid[0] = PHOS Module number 1:fNModules
185 // relid[1] = 0 if PbW04
710f859a 186 // = -1 if CPV
187 // relid[2] = Row number inside a PHOS module
188 // relid[3] = Column number inside a PHOS module
d15a28e7 189
190 Bool_t rv = kTRUE ;
92862013 191 Float_t id = AbsId ;
d15a28e7 192
710f859a 193 Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / GetNCristalsInModule() ) ;
d15a28e7 194
710f859a 195 if ( phosmodulenumber > GetNModules() ) { // it is a CPV pad
196
197 id -= GetNPhi() * GetNZ() * GetNModules() ;
198 Float_t nCPV = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ;
199 relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ;
200 relid[1] = -1 ;
201 id -= ( relid[0] - 1 ) * nCPV ;
202 relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ;
203 relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ;
d15a28e7 204 }
710f859a 205 else { // it is a PW04 crystal
d15a28e7 206
92862013 207 relid[0] = phosmodulenumber ;
208 relid[1] = 0 ;
209 id -= ( phosmodulenumber - 1 ) * GetNPhi() * GetNZ() ;
710f859a 210 relid[2] = (Int_t)TMath::Ceil( id / GetNZ() ) ;
211 relid[3] = (Int_t)( id - ( relid[2] - 1 ) * GetNZ() ) ;
d15a28e7 212 }
213 return rv ;
214}
52a36ffd 215
9f616d61 216//____________________________________________________________________________
7b7c1533 217void AliPHOSGeometry::EmcModuleCoverage(const Int_t mod, Double_t & tm, Double_t & tM, Double_t & pm, Double_t & pM, Option_t * opt) const
9f616d61 218{
a4e98857 219 // calculates the angular coverage in theta and phi of one EMC (=PHOS) module
9f616d61 220
221 Double_t conv ;
cf0c2bc1 222 if ( opt == Radian() )
9f616d61 223 conv = 1. ;
cf0c2bc1 224 else if ( opt == Degre() )
9f616d61 225 conv = 180. / TMath::Pi() ;
226 else {
21cd0c07 227 Warning("EmcModuleCoverage", "%s unknown option; result in radian", opt) ;
9f616d61 228 conv = 1. ;
229 }
230
710f859a 231 Float_t phi = GetPHOSAngle(mod) * (TMath::Pi() / 180.) ;
232 Float_t y0 = GetIPtoCrystalSurface() ;
5868a791 233 Float_t x0 = GetCellStep()*GetNPhi() ;
234 Float_t z0 = GetCellStep()*GetNZ();
235 Double_t angle = TMath::ATan( x0 / y0 / 2 ) ;
710f859a 236 phi = phi + 1.5 * TMath::Pi() ; // to follow the convention of the particle generator(PHOS is between 220 and 320 deg.)
92862013 237 Double_t max = phi - angle ;
238 Double_t min = phi + angle ;
239 pM = TMath::Max(max, min) * conv ;
240 pm = TMath::Min(max, min) * conv ;
9f616d61 241
5868a791 242 angle = TMath::ATan( z0 / y0 / 2 ) ;
92862013 243 max = TMath::Pi() / 2. + angle ; // to follow the convention of the particle generator(PHOS is at 90 deg.)
244 min = TMath::Pi() / 2. - angle ;
245 tM = TMath::Max(max, min) * conv ;
246 tm = TMath::Min(max, min) * conv ;
9f616d61 247
248}
249
250//____________________________________________________________________________
7b7c1533 251void AliPHOSGeometry::EmcXtalCoverage(Double_t & theta, Double_t & phi, Option_t * opt) const
9f616d61 252{
a4e98857 253 // calculates the angular coverage in theta and phi of a single crystal in a EMC(=PHOS) module
9f616d61 254
255 Double_t conv ;
cf0c2bc1 256 if ( opt == Radian() )
9f616d61 257 conv = 1. ;
cf0c2bc1 258 else if ( opt == Degre() )
9f616d61 259 conv = 180. / TMath::Pi() ;
260 else {
21cd0c07 261 Warning("EmcXtalCoverage", "%s unknown option; result in radian", opt) ;
9f616d61 262 conv = 1. ;
263 }
264
710f859a 265 Float_t y0 = GetIPtoCrystalSurface() ;
92862013 266 theta = 2 * TMath::ATan( GetCrystalSize(2) / (2 * y0) ) * conv ;
267 phi = 2 * TMath::ATan( GetCrystalSize(0) / (2 * y0) ) * conv ;
9f616d61 268}
269
270
271//____________________________________________________________________________
52a36ffd 272void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TMatrix & gmat) const
d15a28e7 273{
a4e98857 274 // Calculates the coordinates of a RecPoint and the error matrix in the ALICE global coordinate system
b2a60966 275
d15a28e7 276 AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;
92862013 277 TVector3 localposition ;
d15a28e7 278
279 tmpPHOS->GetLocalPosition(gpos) ;
280
281
282 if ( tmpPHOS->IsEmc() ) // it is a EMC crystal
710f859a 283 { gpos.SetY( - GetIPtoCrystalSurface()) ;
d15a28e7 284
285 }
286 else
710f859a 287 { // it is a CPV
288 gpos.SetY(- GetIPtoUpperCPVsurface() ) ;
d15a28e7 289 }
290
92862013 291 Float_t phi = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ;
292 Double_t const kRADDEG = 180.0 / kPI ;
293 Float_t rphi = phi / kRADDEG ;
d15a28e7 294
92862013 295 TRotation rot ;
296 rot.RotateZ(-rphi) ; // a rotation around Z by angle
d15a28e7 297
92862013 298 TRotation dummy = rot.Invert() ; // to transform from original frame to rotate frame
299 gpos.Transform(rot) ; // rotate the baby
6ad0bfa0 300
d15a28e7 301}
302
303//____________________________________________________________________________
5cda30f6 304void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos) const
d15a28e7 305{
a4e98857 306 // Calculates the coordinates of a RecPoint in the ALICE global coordinate system
b2a60966 307
d15a28e7 308 AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;
92862013 309 TVector3 localposition ;
d15a28e7 310 tmpPHOS->GetLocalPosition(gpos) ;
311
312
313 if ( tmpPHOS->IsEmc() ) // it is a EMC crystal
710f859a 314 { gpos.SetY( - GetIPtoCrystalSurface() ) ;
d15a28e7 315 }
316 else
710f859a 317 { // it is a CPV
318 gpos.SetY(- GetIPtoUpperCPVsurface() ) ;
d15a28e7 319 }
320
92862013 321 Float_t phi = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ;
322 Double_t const kRADDEG = 180.0 / kPI ;
323 Float_t rphi = phi / kRADDEG ;
d15a28e7 324
92862013 325 TRotation rot ;
326 rot.RotateZ(-rphi) ; // a rotation around Z by angle
d15a28e7 327
92862013 328 TRotation dummy = rot.Invert() ; // to transform from original frame to rotate frame
329 gpos.Transform(rot) ; // rotate the baby
d15a28e7 330}
331
332//____________________________________________________________________________
7b7c1533 333void AliPHOSGeometry::ImpactOnEmc(const Double_t theta, const Double_t phi, Int_t & ModuleNumber, Double_t & z, Double_t & x) const
d15a28e7 334{
a4e98857 335 // calculates the impact coordinates on PHOS of a neutral particle
336 // emitted in the direction theta and phi in the ALICE global coordinate system
d15a28e7 337
382f5715 338 //Convert phi to range 0-2pi if nesassary
339 Double_t phiin2pi = phi ;
340 while(phiin2pi<0)
341 phiin2pi+=6.2831853072 ;
342
52a36ffd 343 // searches for the PHOS EMC module
344 ModuleNumber = 0 ;
345 Double_t tm, tM, pm, pM ;
346 Int_t index = 1 ;
347 while ( ModuleNumber == 0 && index <= GetNModules() ) {
348 EmcModuleCoverage(index, tm, tM, pm, pM) ;
382f5715 349 if ( (theta >= tm && theta <= tM) && (phiin2pi >= pm && phiin2pi <= pM ) )
52a36ffd 350 ModuleNumber = index ;
351 index++ ;
d15a28e7 352 }
52a36ffd 353 if ( ModuleNumber != 0 ) {
354 Float_t phi0 = GetPHOSAngle(ModuleNumber) * (TMath::Pi() / 180.) + 1.5 * TMath::Pi() ;
710f859a 355 Float_t y0 = GetIPtoCrystalSurface() ;
382f5715 356 Double_t angle = phiin2pi - phi0;
52a36ffd 357 x = y0 * TMath::Tan(angle) ;
358 angle = theta - TMath::Pi() / 2 ;
359 z = y0 * TMath::Tan(angle) ;
d15a28e7 360 }
d15a28e7 361}
362
0bc3b8ed 363//____________________________________________________________________________
1c9d8212 364Bool_t AliPHOSGeometry::Impact(const TParticle * particle) const
365{
0bc3b8ed 366 // Check if a particle being propagates from IP along the straight line impacts EMC
367
368 Bool_t in=kFALSE;
369 Int_t moduleNumber=0;
1c9d8212 370 Double_t z,x;
0bc3b8ed 371 ImpactOnEmc(particle->Theta(),particle->Phi(),moduleNumber,z,x);
372 if(moduleNumber) in=kTRUE;
373 else in=kFALSE;
374 return in;
1c9d8212 375}
376
d15a28e7 377//____________________________________________________________________________
7b7c1533 378Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t & AbsId) const
d15a28e7 379{
b2a60966 380 // Converts the relative numbering into the absolute numbering
ed4205d8 381 // EMCA crystals:
382 // AbsId = from 1 to fNModules * fNPhi * fNZ
ed4205d8 383 // CPV pad:
384 // AbsId = from N(total PHOS crystals) + 1
385 // to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ
d15a28e7 386
387 Bool_t rv = kTRUE ;
710f859a 388
389 if ( relid[1] == 0 ) { // it is a Phos crystal
52a36ffd 390 AbsId =
710f859a 391 ( relid[0] - 1 ) * GetNPhi() * GetNZ() // the offset of PHOS modules
392 + ( relid[2] - 1 ) * GetNZ() // the offset along phi
393 + relid[3] ; // the offset along z
d15a28e7 394 }
710f859a 395 else { // it is a CPV pad
396 AbsId = GetNPhi() * GetNZ() * GetNModules() // the offset to separate EMCA crystals from CPV pads
397 + ( relid[0] - 1 ) * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() // the pads offset of PHOS modules
398 + ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() // the pads offset of a CPV row
52a36ffd 399 + relid[3] ; // the column number
400 }
401
d15a28e7 402 return rv ;
403}
404
405//____________________________________________________________________________
5868a791 406void AliPHOSGeometry::RelPosToAbsId(const Int_t module , const Double_t x, const Double_t z, Int_t & AbsId)const{
407 // Converts local PHOS-module (x, z) coordinates to absId
d15a28e7 408
5868a791 409 if(!module){
410 AbsId = 0 ;
411 return ;
412 }
413
414 Int_t relid[4] ;
415 relid[0] = module ;
416 relid[1] = 0 ;
417 relid[2] = static_cast<Int_t>(TMath::Ceil(GetNPhi()/2.+ x/GetCellStep()));
418 relid[3] = static_cast<Int_t>(TMath::Ceil(GetNZ()/2. - z/GetCellStep())) ;
419
420 RelToAbsNumbering(relid,AbsId) ;
421
422}
423//____________________________________________________________________________
7b7c1533 424void AliPHOSGeometry::RelPosInAlice(const Int_t id, TVector3 & pos ) const
d15a28e7 425{
a4e98857 426 // Converts the absolute numbering into the global ALICE coordinate system
b2a60966 427
ed4205d8 428
429 Int_t relid[4] ;
430
431 AbsToRelNumbering(id , relid) ;
432
433 Int_t phosmodule = relid[0] ;
434
435 Float_t y0 = 0 ;
436
710f859a 437 if ( relid[1] == 0 ) // it is a PbW04 crystal
438 y0 = - GetIPtoCrystalSurface() ;
439 else
440 y0 = - GetIPtoUpperCPVsurface() ;
441
ed4205d8 442 Float_t x, z ;
443 RelPosInModule(relid, x, z) ;
444
445 pos.SetX(x) ;
446 pos.SetZ(z) ;
710f859a 447 pos.SetY(y0) ;
ed4205d8 448
449 Float_t phi = GetPHOSAngle( phosmodule) ;
450 Double_t const kRADDEG = 180.0 / kPI ;
451 Float_t rphi = phi / kRADDEG ;
452
453 TRotation rot ;
454 rot.RotateZ(-rphi) ; // a rotation around Z by angle
455
456 TRotation dummy = rot.Invert() ; // to transform from original frame to rotate frame
457
458 pos.Transform(rot) ; // rotate the baby
d15a28e7 459}
460
461//____________________________________________________________________________
7b7c1533 462void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) const
d15a28e7 463{
b2a60966 464 // Converts the relative numbering into the local PHOS-module (x, z) coordinates
52a36ffd 465 // Note: sign of z differs from that in the previous version (Yu.Kharlov, 12 Oct 2000)
b2a60966 466
786222b3 467 Int_t row = relid[2] ; //offset along x axis
468 Int_t column = relid[3] ; //offset along z axis
d15a28e7 469
ed4205d8 470
66c3e8ff 471 if ( relid[1] == 0 ) { // its a PbW04 crystal
472 x = - ( GetNPhi()/2. - row + 0.5 ) * GetCellStep() ; // position of Xtal with respect
473 z = ( GetNZ() /2. - column + 0.5 ) * GetCellStep() ; // of center of PHOS module
52a36ffd 474 }
475 else {
710f859a 476 x = - ( GetNumberOfCPVPadsPhi()/2. - row - 0.5 ) * GetPadSizePhi() ; // position of pad with respect
477 z = ( GetNumberOfCPVPadsZ() /2. - column - 0.5 ) * GetPadSizeZ() ; // of center of PHOS module
52a36ffd 478 }
2f3366b6 479}