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