Make separate, specialized geometries for RPhi and RhoZ views.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSGeometry.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 PHOS : singleton  
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.
26 //                  
27 // -- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (RRC "KI" & SUBATECH)
28
29 // --- ROOT system ---
30
31 #include "TVector3.h"
32 #include "TRotation.h" 
33 #include "TParticle.h"
34 #include <TGeoManager.h>
35 #include <TGeoMatrix.h>
36
37 // --- Standard library ---
38
39 // --- AliRoot header files ---
40 #include "AliLog.h"
41 #include "AliPHOSGeometry.h"
42 #include "AliPHOSEMCAGeometry.h" 
43 #include "AliPHOSRecPoint.h"
44
45 ClassImp(AliPHOSGeometry)
46
47 // these initialisations are needed for a singleton
48 AliPHOSGeometry  * AliPHOSGeometry::fgGeom = 0 ;
49 Bool_t             AliPHOSGeometry::fgInit = kFALSE ;
50
51 //____________________________________________________________________________
52 AliPHOSGeometry::AliPHOSGeometry() : 
53                     fNModules(0),
54                     fAngle(0.f),
55                     fPHOSAngle(0),
56                     fIPtoUpperCPVsurface(0),
57                     fCrystalShift(0),
58                     fCryCellShift(0),
59                     fRotMatrixArray(0),
60                     fGeometryEMCA(0),
61                     fGeometryCPV(0),
62                     fGeometrySUPP(0)
63 {
64     // default ctor 
65     // must be kept public for root persistency purposes, but should never be called by the outside world
66     fgGeom          = 0 ;
67 }  
68
69 //____________________________________________________________________________
70 AliPHOSGeometry::AliPHOSGeometry(const AliPHOSGeometry & rhs)
71                     : AliGeometry(rhs),
72                       fNModules(rhs.fNModules),
73                       fAngle(rhs.fAngle),
74                       fPHOSAngle(0),
75                       fIPtoUpperCPVsurface(rhs.fIPtoUpperCPVsurface),
76                       fCrystalShift(rhs.fCrystalShift),
77                       fCryCellShift(rhs.fCryCellShift),
78                       fRotMatrixArray(0),
79                       fGeometryEMCA(0),
80                       fGeometryCPV(0),
81                       fGeometrySUPP(0)
82 {
83   Fatal("cpy ctor", "not implemented") ; 
84 }
85
86 //____________________________________________________________________________
87 AliPHOSGeometry::AliPHOSGeometry(const Text_t* name, const Text_t* title) 
88                   : AliGeometry(name, title),
89                     fNModules(0),
90                     fAngle(0.f),
91                     fPHOSAngle(0),
92                     fIPtoUpperCPVsurface(0),
93                     fCrystalShift(0),
94                     fCryCellShift(0),
95                     fRotMatrixArray(0),
96                     fGeometryEMCA(0),
97                     fGeometryCPV(0),
98                     fGeometrySUPP(0)
99
100   // ctor only for internal usage (singleton)
101   Init() ; 
102   fgGeom = this;
103 }
104
105 //____________________________________________________________________________
106 AliPHOSGeometry::~AliPHOSGeometry(void)
107 {
108   // dtor
109
110   if (fRotMatrixArray) fRotMatrixArray->Delete() ; 
111   if (fRotMatrixArray) delete fRotMatrixArray ; 
112   if (fPHOSAngle     ) delete[] fPHOSAngle ; 
113 }
114
115 //____________________________________________________________________________
116 void AliPHOSGeometry::Init(void)
117 {
118   // Initializes the PHOS parameters :
119   //  IHEP is the Protvino CPV (cathode pad chambers)
120   
121 /*
122   TString test(GetName()) ; 
123   if (test != "IHEP" && test != "noCPV") {
124     AliFatal(Form("%s is not a known geometry (choose among IHEP)", 
125                   test.Data() )) ; 
126   }
127 */
128
129   fgInit     = kTRUE ; 
130
131   fNModules     = 5;
132   fAngle        = 20;
133
134   fGeometryEMCA = new AliPHOSEMCAGeometry();
135   
136   fGeometryCPV  = new AliPHOSCPVGeometry ();
137   
138   fGeometrySUPP = new AliPHOSSupportGeometry();
139   
140   fPHOSAngle = new Float_t[fNModules] ;
141   
142   Float_t * emcParams = fGeometryEMCA->GetEMCParams() ;
143   
144   fPHOSParams[0] =  TMath::Max((Double_t)fGeometryCPV->GetCPVBoxSize(0)/2., 
145                                (Double_t)(emcParams[0] - (emcParams[1]-emcParams[0])*
146                                           fGeometryCPV->GetCPVBoxSize(1)/2/emcParams[3]));
147   fPHOSParams[1] = emcParams[1] ;
148   fPHOSParams[2] = TMath::Max((Double_t)emcParams[2], (Double_t)fGeometryCPV->GetCPVBoxSize(2)/2.);
149   fPHOSParams[3] = emcParams[3] + fGeometryCPV->GetCPVBoxSize(1)/2. ;
150   
151   fIPtoUpperCPVsurface = fGeometryEMCA->GetIPtoOuterCoverDistance() - fGeometryCPV->GetCPVBoxSize(1) ;
152
153   //calculate offset to crystal surface
154   Float_t * inthermo = fGeometryEMCA->GetInnerThermoHalfSize() ;
155   Float_t * strip = fGeometryEMCA->GetStripHalfSize() ;
156   Float_t* splate = fGeometryEMCA->GetSupportPlateHalfSize();
157   Float_t * crystal = fGeometryEMCA->GetCrystalHalfSize() ;
158   Float_t * pin = fGeometryEMCA->GetAPDHalfSize() ;
159   Float_t * preamp = fGeometryEMCA->GetPreampHalfSize() ;
160   fCrystalShift=-inthermo[1]+strip[1]+splate[1]+crystal[1]-fGeometryEMCA->GetAirGapLed()/2.+pin[1]+preamp[1] ;
161   fCryCellShift=crystal[1]-(fGeometryEMCA->GetAirGapLed()-2*pin[1]-2*preamp[1])/2;
162  
163   Int_t index ;
164   for ( index = 0; index < fNModules; index++ )
165     fPHOSAngle[index] = 0.0 ; // Module position angles are set in CreateGeometry()
166   
167   fRotMatrixArray = new TObjArray(fNModules) ; 
168
169   // Geometry parameters are calculated
170
171   SetPHOSAngles();
172   Double_t const kRADDEG = 180.0 / TMath::Pi() ;
173   Float_t r = GetIPtoOuterCoverDistance() + fPHOSParams[3] - GetCPVBoxSize(1) ;
174   for (Int_t iModule=0; iModule<fNModules; iModule++) {
175     fModuleCenter[iModule][0] = r * TMath::Sin(fPHOSAngle[iModule] / kRADDEG );
176     fModuleCenter[iModule][1] =-r * TMath::Cos(fPHOSAngle[iModule] / kRADDEG );
177     fModuleCenter[iModule][2] = 0.;
178     
179     fModuleAngle[iModule][0][0] =  90;
180     fModuleAngle[iModule][0][1] =   fPHOSAngle[iModule];
181     fModuleAngle[iModule][1][0] =   0;
182     fModuleAngle[iModule][1][1] =   0;
183     fModuleAngle[iModule][2][0] =  90;
184     fModuleAngle[iModule][2][1] = 270 + fPHOSAngle[iModule];
185   }
186
187 }
188
189 //____________________________________________________________________________
190 AliPHOSGeometry *  AliPHOSGeometry::GetInstance() 
191
192   // Returns the pointer of the unique instance; singleton specific
193   
194   return static_cast<AliPHOSGeometry *>( fgGeom ) ; 
195 }
196
197 //____________________________________________________________________________
198 AliPHOSGeometry *  AliPHOSGeometry::GetInstance(const Text_t* name, const Text_t* title) 
199 {
200   // Returns the pointer of the unique instance
201   // Creates it with the specified options (name, title) if it does not exist yet
202
203   AliPHOSGeometry * rv = 0  ; 
204   if ( fgGeom == 0 ) {
205     if ( strcmp(name,"") == 0 ) 
206       rv = 0 ;
207     else {    
208       fgGeom = new AliPHOSGeometry(name, title) ;
209       if ( fgInit )
210         rv = (AliPHOSGeometry * ) fgGeom ;
211       else {
212         rv = 0 ; 
213         delete fgGeom ; 
214         fgGeom = 0 ; 
215       }
216     }
217   }
218   else {
219     if ( strcmp(fgGeom->GetName(), name) != 0 ) 
220       ::Error("GetInstance", "Current geometry is %s. You cannot call %s", 
221                       fgGeom->GetName(), name) ; 
222     else
223       rv = (AliPHOSGeometry *) fgGeom ; 
224   } 
225   return rv ; 
226 }
227
228 //____________________________________________________________________________
229 void AliPHOSGeometry::SetPHOSAngles() 
230
231   // Calculates the position of the PHOS modules in ALICE global coordinate system
232   // in ideal geometry
233   
234   Double_t const kRADDEG = 180.0 / TMath::Pi() ;
235   Float_t pphi =  2 * TMath::ATan( GetOuterBoxSize(0)  / ( 2.0 * GetIPtoUpperCPVsurface() ) ) ;
236   pphi *= kRADDEG ;
237   if (pphi > fAngle){ 
238     AliError(Form("PHOS modules overlap!\n pphi = %f fAngle = %f", 
239                   pphi, fAngle));
240
241   }
242   pphi = fAngle;
243   
244   for( Int_t i = 1; i <= fNModules ; i++ ) {
245     Float_t angle = pphi * ( i - fNModules / 2.0 - 0.5 ) ;
246     fPHOSAngle[i-1] = -  angle ;
247   } 
248 }
249
250 //____________________________________________________________________________
251 Bool_t AliPHOSGeometry::AbsToRelNumbering(Int_t absId, Int_t * relid) const
252 {
253   // Converts the absolute numbering into the following array
254   //  relid[0] = PHOS Module number 1:fNModules 
255   //  relid[1] = 0 if PbW04
256   //           = -1 if CPV
257   //  relid[2] = Row number inside a PHOS module
258   //  relid[3] = Column number inside a PHOS module
259
260   Bool_t rv  = kTRUE ; 
261   Float_t id = absId ;
262
263   Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / GetNCristalsInModule() ) ; 
264   
265   if ( phosmodulenumber >  GetNModules() ) { // it is a CPV pad
266     
267     id -=  GetNPhi() * GetNZ() *  GetNModules() ; 
268     Float_t nCPV  = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ;
269     relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ;
270     relid[1] = -1 ;
271     id -= ( relid[0] - 1 ) * nCPV ; 
272     relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ;
273     relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ; 
274   } 
275   else { // it is a PW04 crystal
276
277     relid[0] = phosmodulenumber ;
278     relid[1] = 0 ;
279     id -= ( phosmodulenumber - 1 ) *  GetNPhi() * GetNZ() ; 
280     relid[2] = (Int_t)TMath::Ceil( id / GetNZ() )  ;
281     relid[3] = (Int_t)( id - ( relid[2] - 1 ) * GetNZ() ) ; 
282   } 
283   return rv ; 
284 }
285 //____________________________________________________________________________
286 void AliPHOSGeometry::GetGlobal(const AliRecPoint* , TVector3 & ) const
287 {
288   AliFatal(Form("Please use GetGlobalPHOS(recPoint,gpos) instead of GetGlobal!"));
289 }
290
291 //____________________________________________________________________________
292 void AliPHOSGeometry::GetGlobalPHOS(const AliPHOSRecPoint* recPoint, TVector3 & gpos) const
293 {
294   // Calculates the coordinates of a RecPoint and the error matrix in the ALICE global coordinate system
295  
296   const AliPHOSRecPoint * tmpPHOS = recPoint ;  
297   TVector3 localposition ;
298
299   tmpPHOS->GetLocalPosition(gpos) ;
300
301   if (!gGeoManager){
302     AliFatal("Geo manager not initialized\n");
303   }
304   //construct module name
305   char path[100] ; 
306   Double_t dy ;
307   if(tmpPHOS->IsEmc()){
308     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",tmpPHOS->GetPHOSMod()) ;
309 //    sprintf(path,"/ALIC_1/PHOS_%d",tmpPHOS->GetPHOSMod()) ;
310     dy=fCrystalShift ;
311   }
312   else{
313     sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",tmpPHOS->GetPHOSMod()) ;
314     dy= GetCPVBoxSize(1)/2. ; //center of CPV module 
315   }
316   Double_t pos[3]={gpos.X(),gpos.Y()-dy,gpos.Z()} ;
317   if(tmpPHOS->IsEmc())
318     pos[2]=-pos[2] ; //Opposite z directions in EMC matrix and local frame!!!
319   Double_t posC[3];
320   //now apply possible shifts and rotations
321   if (!gGeoManager->cd(path)){
322     AliFatal("Geo manager can not find path \n");
323   }
324   TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
325   if (m){
326      m->LocalToMaster(pos,posC);
327   }
328   else{
329     AliFatal("Geo matrixes are not loaded \n") ;
330   }
331   gpos.SetXYZ(posC[0],posC[1],posC[2]) ;
332
333 }
334 //____________________________________________________________________________
335 void AliPHOSGeometry::ImpactOnEmc(Double_t * vtx, Double_t theta, Double_t phi, 
336                                   Int_t & moduleNumber, Double_t & z, Double_t & x) const
337 {
338   // calculates the impact coordinates on PHOS of a neutral particle  
339   // emitted in the vertex vtx[3] with direction theta and phi in the ALICE global coordinate system
340   TVector3 p(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) ;
341   TVector3 v(vtx[0],vtx[1],vtx[2]) ;
342
343   if (!gGeoManager){
344     AliFatal("Geo manager not initialized\n");
345   }
346  
347   for(Int_t imod=1; imod<=GetNModules() ; imod++){
348     //create vector from (0,0,0) to center of crystal surface of imod module
349     Double_t tmp[3]={0.,-fCrystalShift,0.} ;
350  
351     char path[100] ;
352     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",imod) ;
353     if (!gGeoManager->cd(path)){
354       AliFatal("Geo manager can not find path \n");
355     }
356     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
357     Double_t posG[3]={0.,0.,0.} ;
358     if (m) m->LocalToMaster(tmp,posG);
359     TVector3 n(posG[0],posG[1],posG[2]) ; 
360     Double_t direction=n.Dot(p) ;
361     if(direction<=0.)
362       continue ; //momentum directed FROM module
363     Double_t fr = (n.Mag2()-n.Dot(v))/direction ;  
364     //Calculate direction in module plain
365     n-=v+fr*p ;
366     n*=-1. ;
367     Float_t * sz = fGeometryEMCA->GetInnerThermoHalfSize() ; //Wery close to the zise of the Xtl set
368     if(TMath::Abs(TMath::Abs(n.Z())<sz[2]) && n.Pt()<sz[0]){
369       moduleNumber = imod ;
370       z=n.Z() ;
371       x=TMath::Sign(n.Pt(),n.X()) ;
372       //no need to return to local system since we calcilated distance from module center
373       //and tilts can not be significant.
374       return ;
375     }
376   }
377   //Not in acceptance
378   x=0; z=0 ;
379   moduleNumber=0 ;
380
381 }
382
383 //____________________________________________________________________________
384 Bool_t  AliPHOSGeometry::Impact(const TParticle * particle) const 
385 {
386   // Tells if a particle enters PHOS
387   Bool_t in=kFALSE;
388   Int_t moduleNumber=0;
389   Double_t vtx[3]={particle->Vx(),particle->Vy(),particle->Vz()} ;
390   Double_t z,x;
391   ImpactOnEmc(vtx,particle->Theta(),particle->Phi(),moduleNumber,z,x);
392   if(moduleNumber!=0) 
393     in=kTRUE;
394   return in;
395 }
396
397 //____________________________________________________________________________
398 Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  absId) const
399 {
400   // Converts the relative numbering into the absolute numbering
401   // EMCA crystals:
402   //  absId = from 1 to fNModules * fNPhi * fNZ
403   // CPV pad:
404   //  absId = from N(total PHOS crystals) + 1
405   //          to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ
406
407   Bool_t rv = kTRUE ; 
408   
409   if ( relid[1] ==  0 ) {                            // it is a Phos crystal
410     absId =
411       ( relid[0] - 1 ) * GetNPhi() * GetNZ()         // the offset of PHOS modules
412       + ( relid[2] - 1 ) * GetNZ()                   // the offset along phi
413       +   relid[3] ;                                 // the offset along z
414   }
415   else { // it is a CPV pad
416     absId =    GetNPhi() * GetNZ() *  GetNModules()         // the offset to separate EMCA crystals from CPV pads
417       + ( relid[0] - 1 ) * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ()   // the pads offset of PHOS modules 
418       + ( relid[2] - 1 ) * GetNumberOfCPVPadsZ()                             // the pads offset of a CPV row
419       +   relid[3] ;                                                         // the column number
420   }
421   
422   return rv ; 
423 }
424
425 //____________________________________________________________________________
426 void AliPHOSGeometry::RelPosInAlice(Int_t id, TVector3 & pos ) const
427 {
428   // Converts the absolute numbering into the global ALICE coordinate system
429   
430   if (!gGeoManager){
431     AliFatal("Geo manager not initialized\n");
432   }
433     
434   Int_t relid[4] ;
435     
436   AbsToRelNumbering(id , relid) ;
437     
438   //construct module name
439   char path[100] ;
440   if(relid[1]==0){ //this is EMC
441  
442     Double_t ps[3]= {0.0,-fCryCellShift,0.}; //Position incide the crystal 
443     Double_t psC[3]={0.0,0.0,0.}; //Global position
444  
445     //Shift and possibly apply misalignment corrections
446     Int_t nCellsXInStrip=fGeometryEMCA->GetNCellsXInStrip() ;
447     Int_t nCellsZInStrip=fGeometryEMCA->GetNCellsZInStrip() ;
448     Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/nCellsXInStrip))*fGeometryEMCA->GetNStripZ()-
449                 (Int_t) TMath::Ceil((Double_t)relid[3]/nCellsZInStrip) ;
450     Int_t cellraw= relid[3]%nCellsZInStrip ;
451     if(cellraw==0)cellraw=nCellsZInStrip ;
452     Int_t cell= ((relid[2]-1)%nCellsXInStrip)*nCellsZInStrip + cellraw ;
453     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1/PSTR_%d/PCEL_%d",
454             relid[0],strip,cell) ;
455     if (!gGeoManager->cd(path)){
456       AliFatal("Geo manager can not find path \n");
457     }
458     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
459     if (m) m->LocalToMaster(ps,psC);
460     else{
461       AliFatal("Geo matrixes are not loaded \n") ;
462     }
463     pos.SetXYZ(psC[0],psC[1],psC[2]) ; 
464   }
465   else{
466     //first calculate position with respect to CPV plain
467     Int_t row        = relid[2] ; //offset along x axis
468     Int_t column     = relid[3] ; //offset along z axis
469     Double_t ps[3]= {0.0,GetCPVBoxSize(1)/2.,0.}; //Position on top of CPV
470     Double_t psC[3]={0.0,0.0,0.}; //Global position
471     pos[0] = - ( GetNumberOfCPVPadsPhi()/2. - row    - 0.5 ) * GetPadSizePhi()  ; // position of pad  with respect
472     pos[2] = - ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module
473  
474     //now apply possible shifts and rotations
475     sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",relid[0]) ;
476     if (!gGeoManager->cd(path)){
477       AliFatal("Geo manager can not find path \n");
478     }
479     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
480     if (m) m->LocalToMaster(ps,psC);
481     else{
482       AliFatal("Geo matrixes are not loaded \n") ;
483     }
484     pos.SetXYZ(psC[0],psC[1],-psC[2]) ; 
485   }
486
487
488 //____________________________________________________________________________
489 void AliPHOSGeometry::RelPosToAbsId(Int_t module, Double_t x, Double_t z, Int_t & absId) const
490 {
491   // converts local PHOS-module (x, z) coordinates to absId 
492
493   //find Global position
494   if (!gGeoManager){
495     AliFatal("Geo manager not initialized\n");
496   }
497   Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!!
498   Double_t posG[3] ;
499   char path[100] ;
500   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",module) ;
501   if (!gGeoManager->cd(path)){
502     AliFatal("Geo manager can not find path \n");
503   }
504   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
505   if (mPHOS){
506      mPHOS->LocalToMaster(posL,posG);
507   }
508   else{
509     AliFatal("Geo matrixes are not loaded \n") ;
510   }
511
512   Int_t relid[4] ;
513   gGeoManager->FindNode(posG[0],posG[1],posG[2]) ;
514   //Check that path contains PSTR and extract strip number
515   TString cpath(gGeoManager->GetPath()) ;
516   Int_t indx = cpath.Index("PCEL") ;
517   if(indx==-1){ //for the few events when particle hits between srips use ideal geometry
518     relid[0] = module ;
519     relid[1] = 0 ;
520     relid[2] = static_cast<Int_t>(TMath::Ceil( x/ GetCellStep() + GetNPhi() / 2.) );
521     relid[3] = static_cast<Int_t>(TMath::Ceil(-z/ GetCellStep() + GetNZ()   / 2.) ) ;
522     if(relid[2]<1)relid[2]=1 ;
523     if(relid[3]<1)relid[3]=1 ;
524     if(relid[2]>GetNPhi())relid[2]=GetNPhi() ;
525     if(relid[3]>GetNZ())relid[3]=GetNZ() ;
526     RelToAbsNumbering(relid,absId) ;
527   }
528   else{
529     Int_t indx2 = cpath.Index("/",indx) ;
530     if(indx2==-1)
531       indx2=cpath.Length() ;
532     TString cell=cpath(indx+5,indx2-indx-5) ;
533     Int_t icell=cell.Atoi() ;
534     indx = cpath.Index("PSTR") ;
535     indx2 = cpath.Index("/",indx) ;
536     TString strip=cpath(indx+5,indx2-indx-5) ;
537     Int_t iStrip = strip.Atoi() ; 
538
539     Int_t row = fGeometryEMCA->GetNStripZ() - (iStrip - 1) % (fGeometryEMCA->GetNStripZ()) ;
540     Int_t col = (Int_t) TMath::Ceil((Double_t) iStrip/(fGeometryEMCA->GetNStripZ())) -1 ;
541  
542     // Absid for 8x2-strips. Looks nice :)
543     absId = (module-1)*GetNCristalsInModule() +
544                   row * 2 + (col*fGeometryEMCA->GetNCellsXInStrip() + (icell - 1) / 2)*GetNZ() - (icell & 1 ? 1 : 0);
545  
546   }
547  
548 }
549
550 //____________________________________________________________________________
551 void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) const 
552 {
553   // Converts the relative numbering into the local PHOS-module (x, z) coordinates
554   // Note: sign of z differs from that in the previous version (Yu.Kharlov, 12 Oct 2000)
555   
556
557   if (!gGeoManager){
558     AliFatal("Geo manager not initialized\n");
559   }
560   //construct module name
561   char path[100] ;
562   if(relid[1]==0){ //this is PHOS
563
564 //   Calculations using ideal geometry (obsolete)
565 //    x = - ( GetNPhi()/2. - relid[2]    + 0.5 ) *  GetCellStep() ; // position of Xtal with respect
566 //    z = - ( GetNZ()  /2. - relid[3] + 0.5 ) *  GetCellStep() ; // of center of PHOS module  
567
568     Double_t pos[3]= {0.0,-fCryCellShift,0.}; //Position incide the crystal 
569     Double_t posC[3]={0.0,0.0,0.}; //Global position
570
571     //Shift and possibly apply misalignment corrections
572     Int_t nCellsXInStrip=fGeometryEMCA->GetNCellsXInStrip() ;
573     Int_t nCellsZInStrip=fGeometryEMCA->GetNCellsZInStrip() ;
574 //    Int_t strip=1+(relid[3]-1)/fGeometryEMCA->GetNCellsZInStrip()+((relid[2]-1)/nCellsInStrip)*fGeometryEMCA->GetNStripZ() ;
575     Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/nCellsXInStrip))*fGeometryEMCA->GetNStripZ()-
576                 (Int_t) TMath::Ceil((Double_t)relid[3]/nCellsZInStrip) ;
577     Int_t cellraw= relid[3]%nCellsZInStrip ;
578     if(cellraw==0)cellraw=nCellsZInStrip ;
579     Int_t cell= ((relid[2]-1)%nCellsXInStrip)*nCellsZInStrip + cellraw ; 
580     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1/PSTR_%d/PCEL_%d",
581             relid[0],strip,cell) ;
582     if (!gGeoManager->cd(path)){
583       AliFatal("Geo manager can not find path \n");
584     }
585     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
586     if (m) m->LocalToMaster(pos,posC);
587     else{
588       AliFatal("Geo matrixes are not loaded \n") ;
589     }
590     //    printf("Local: x=%f, y=%f, z=%f \n",pos[0],pos[1],pos[2]) ;
591     //    printf("   gl: x=%f, y=%f, z=%f \n",posC[0],posC[1],posC[2]) ;
592     //Return to PHOS local system  
593     Double_t posL[3]={posC[0],posC[1],posC[2]};
594     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",relid[0]) ;
595     //    sprintf(path,"/ALIC_1/PHOS_%d",relid[0]) ;
596     if (!gGeoManager->cd(path)){
597       AliFatal("Geo manager can not find path \n");
598     }
599     TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
600     if (mPHOS) mPHOS->MasterToLocal(posC,posL);
601     else{
602       AliFatal("Geo matrixes are not loaded \n") ;
603     }
604 //printf("RelPosInMod: posL=[%f,%f,%f]\n",posL[0],posL[1],posL[2]) ;
605 //printf("old: x=%f, z=%f \n",x,z);
606     x=posL[0] ;
607     z=-posL[2];
608     return ;
609   }
610   else{//CPV
611     //first calculate position with respect to CPV plain 
612     Int_t row        = relid[2] ; //offset along x axis
613     Int_t column     = relid[3] ; //offset along z axis
614     Double_t pos[3]= {0.0,0.0,0.}; //Position incide the CPV printed circuit
615     Double_t posC[3]={0.0,0.0,0.}; //Global position
616     //    x = - ( GetNumberOfCPVPadsPhi()/2. - row    - 0.5 ) * GetPadSizePhi()  ; // position of pad  with respect
617     //    z = - ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module
618     pos[0] = - ( GetNumberOfCPVPadsPhi()/2. - row    - 0.5 ) * GetPadSizePhi()  ; // position of pad  with respect
619     pos[2] = - ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module
620
621     //now apply possible shifts and rotations
622     sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",relid[0]) ;
623     if (!gGeoManager->cd(path)){
624       AliFatal("Geo manager can not find path \n");
625     }
626     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
627     if (m) m->LocalToMaster(pos,posC);
628     else{
629       AliFatal("Geo matrixes are not loaded \n") ;
630     }
631     //Return to PHOS local system
632     Double_t posL[3]={0.,0.,0.,} ;
633     sprintf(path,"/ALIC_1/PHOS_%d",relid[0]) ;
634     if (!gGeoManager->cd(path)){
635       AliFatal("Geo manager can not find path \n");
636     }
637     TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
638     if (mPHOS) mPHOS->MasterToLocal(posC,posL);
639     else{
640       AliFatal("Geo matrixes are not loaded \n") ;
641     }
642     x=posL[0] ;
643     z=posL[1];
644     return ;
645  
646   }
647   
648 }
649
650 //____________________________________________________________________________
651
652 void AliPHOSGeometry::GetModuleCenter(TVector3& center, 
653                                       const char *det,
654                                       Int_t module) const
655 {
656   // Returns a position of the center of the CPV or EMC module
657   // in ideal (not misaligned) geometry
658   Float_t rDet = 0.;
659   if      (strcmp(det,"CPV") == 0) rDet  = GetIPtoCPVDistance   ();
660   else if (strcmp(det,"EMC") == 0) rDet  = GetIPtoCrystalSurface();
661   else 
662     AliFatal(Form("Wrong detector name %s",det));
663
664   Float_t angle = GetPHOSAngle(module); // (40,20,0,-20,-40) degrees
665   angle *= TMath::Pi()/180;
666   angle += 3*TMath::Pi()/2.;
667   center.SetXYZ(rDet*TMath::Cos(angle), rDet*TMath::Sin(angle), 0.);
668 }
669
670 //____________________________________________________________________________
671
672 void AliPHOSGeometry::Global2Local(TVector3& localPosition,
673                                    const TVector3& globalPosition,
674                                    Int_t module) const
675 {
676   // Transforms a global position of the rec.point to the local coordinate system
677   //Return to PHOS local system
678   Double_t posG[3]={globalPosition.X(),globalPosition.Y(),globalPosition.Z()} ;
679   Double_t posL[3]={0.,0.,0.} ;
680   char path[100] ;
681   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",module) ;
682 //  sprintf(path,"/ALIC_1/PHOS_%d",module) ;
683   if (!gGeoManager->cd(path)){
684     AliFatal("Geo manager can not find path \n");
685   }
686   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
687   if (mPHOS) mPHOS->MasterToLocal(posG,posL);
688   else{
689     AliFatal("Geo matrixes are not loaded \n") ;
690   }
691   localPosition.SetXYZ(posL[0],posL[1]+fCrystalShift,-posL[2]) ;  
692  
693 /*
694   Float_t angle = GetPHOSAngle(module); // (40,20,0,-20,-40) degrees
695   angle *= TMath::Pi()/180;
696   angle += 3*TMath::Pi()/2.;
697   localPosition = globalPosition;
698   localPosition.RotateZ(-angle);
699 */
700 }
701 //____________________________________________________________________________
702 void AliPHOSGeometry::Local2Global(Int_t mod, Float_t x, Float_t z,
703                                    TVector3& globalPosition) const 
704 {
705   char path[100] ;
706   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",mod) ;
707 //  sprintf(path,"/ALIC_1/PHOS_%d",mod) ;
708   if (!gGeoManager->cd(path)){
709     AliFatal("Geo manager can not find path \n");
710   }
711   Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!!
712   Double_t posG[3] ;
713   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
714   if (mPHOS){
715      mPHOS->LocalToMaster(posL,posG);
716   }    
717   else{
718     AliFatal("Geo matrixes are not loaded \n") ;
719   }
720   globalPosition.SetXYZ(posG[0],posG[1],posG[2]) ;
721 }
722 //____________________________________________________________________________
723 void AliPHOSGeometry::GetIncidentVector(const TVector3 &vtx, Int_t module, Float_t x,Float_t z, TVector3 &vInc) const {
724   //Calculates vector pointing from vertex to current poisition in module local frame
725   //Note that PHOS local system and ALICE global have opposite z directions
726
727   Global2Local(vInc,vtx,module) ; 
728   vInc.SetXYZ(vInc.X()+x,vInc.Y(),vInc.Z()+z) ;
729 }