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