Loaders removed from the reconstruction code (C.Cheshkov)
[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       return ;
364     }
365   }
366   //Not in acceptance
367   x=0; z=0 ;
368   moduleNumber=0 ;
369
370 }
371
372 //____________________________________________________________________________
373 Bool_t  AliPHOSGeometry::Impact(const TParticle * particle) const 
374 {
375   // Tells if a particle enters PHOS
376   Bool_t in=kFALSE;
377   Int_t moduleNumber=0;
378   Double_t vtx[3]={particle->Vx(),particle->Vy(),particle->Vz()} ;
379   Double_t z,x;
380   ImpactOnEmc(vtx,particle->Theta(),particle->Phi(),moduleNumber,z,x);
381   if(moduleNumber!=0) 
382     in=kTRUE;
383   return in;
384 }
385
386 //____________________________________________________________________________
387 Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  absId) const
388 {
389   // Converts the relative numbering into the absolute numbering
390   // EMCA crystals:
391   //  absId = from 1 to fNModules * fNPhi * fNZ
392   // CPV pad:
393   //  absId = from N(total PHOS crystals) + 1
394   //          to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ
395
396   Bool_t rv = kTRUE ; 
397   
398   if ( relid[1] ==  0 ) {                            // it is a Phos crystal
399     absId =
400       ( relid[0] - 1 ) * GetNPhi() * GetNZ()         // the offset of PHOS modules
401       + ( relid[2] - 1 ) * GetNZ()                   // the offset along phi
402       +   relid[3] ;                                 // the offset along z
403   }
404   else { // it is a CPV pad
405     absId =    GetNPhi() * GetNZ() *  GetNModules()         // the offset to separate EMCA crystals from CPV pads
406       + ( relid[0] - 1 ) * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ()   // the pads offset of PHOS modules 
407       + ( relid[2] - 1 ) * GetNumberOfCPVPadsZ()                             // the pads offset of a CPV row
408       +   relid[3] ;                                                         // the column number
409   }
410   
411   return rv ; 
412 }
413
414 //____________________________________________________________________________
415 void AliPHOSGeometry::RelPosInAlice(Int_t id, TVector3 & pos ) const
416 {
417   // Converts the absolute numbering into the global ALICE coordinate system
418   
419   if (!gGeoManager){
420     AliFatal("Geo manager not initialized\n");
421   }
422     
423   Int_t relid[4] ;
424     
425   AbsToRelNumbering(id , relid) ;
426     
427   //construct module name
428   char path[100] ;
429   if(relid[1]==0){ //this is EMC
430  
431     Double_t ps[3]= {0.0,-fCryCellShift,0.}; //Position incide the crystal 
432     Double_t psC[3]={0.0,0.0,0.}; //Global position
433  
434     //Shift and possibly apply misalignment corrections
435     Int_t nCellsXInStrip=fGeometryEMCA->GetNCellsXInStrip() ;
436     Int_t nCellsZInStrip=fGeometryEMCA->GetNCellsZInStrip() ;
437     Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/nCellsXInStrip))*fGeometryEMCA->GetNStripZ()-
438                 (Int_t) TMath::Ceil((Double_t)relid[3]/nCellsZInStrip) ;
439     Int_t cellraw= relid[3]%nCellsZInStrip ;
440     if(cellraw==0)cellraw=nCellsZInStrip ;
441     Int_t cell= ((relid[2]-1)%nCellsXInStrip)*nCellsZInStrip + cellraw ;
442     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1/PSTR_%d/PCEL_%d",
443             relid[0],strip,cell) ;
444     if (!gGeoManager->cd(path)){
445       AliFatal("Geo manager can not find path \n");
446     }
447     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
448     if (m) m->LocalToMaster(ps,psC);
449     else{
450       AliFatal("Geo matrixes are not loaded \n") ;
451     }
452     pos.SetXYZ(psC[0],psC[1],psC[2]) ; 
453   }
454   else{
455     //first calculate position with respect to CPV plain
456     Int_t row        = relid[2] ; //offset along x axis
457     Int_t column     = relid[3] ; //offset along z axis
458     Double_t ps[3]= {0.0,GetCPVBoxSize(1)/2.,0.}; //Position on top of CPV
459     Double_t psC[3]={0.0,0.0,0.}; //Global position
460     pos[0] = - ( GetNumberOfCPVPadsPhi()/2. - row    - 0.5 ) * GetPadSizePhi()  ; // position of pad  with respect
461     pos[2] = - ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module
462  
463     //now apply possible shifts and rotations
464     sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",relid[0]) ;
465     if (!gGeoManager->cd(path)){
466       AliFatal("Geo manager can not find path \n");
467     }
468     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
469     if (m) m->LocalToMaster(ps,psC);
470     else{
471       AliFatal("Geo matrixes are not loaded \n") ;
472     }
473     pos.SetXYZ(psC[0],psC[1],-psC[2]) ; 
474   }
475
476
477 //____________________________________________________________________________
478 void AliPHOSGeometry::RelPosToAbsId(Int_t module, Double_t x, Double_t z, Int_t & absId) const
479 {
480   // converts local PHOS-module (x, z) coordinates to absId 
481
482   //find Global position
483   if (!gGeoManager){
484     AliFatal("Geo manager not initialized\n");
485   }
486   Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!!
487   Double_t posG[3] ;
488   char path[100] ;
489   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",module) ;
490   if (!gGeoManager->cd(path)){
491     AliFatal("Geo manager can not find path \n");
492   }
493   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
494   if (mPHOS){
495      mPHOS->LocalToMaster(posL,posG);
496   }
497   else{
498     AliFatal("Geo matrixes are not loaded \n") ;
499   }
500
501   Int_t relid[4] ;
502   gGeoManager->FindNode(posG[0],posG[1],posG[2]) ;
503   //Check that path contains PSTR and extract strip number
504   TString cpath(gGeoManager->GetPath()) ;
505   Int_t indx = cpath.Index("PCEL") ;
506   if(indx==-1){ //for the few events when particle hits between srips use ideal geometry
507     relid[0] = module ;
508     relid[1] = 0 ;
509     relid[2] = static_cast<Int_t>(TMath::Ceil( x/ GetCellStep() + GetNPhi() / 2.) );
510     relid[3] = static_cast<Int_t>(TMath::Ceil(-z/ GetCellStep() + GetNZ()   / 2.) ) ;
511     if(relid[2]<1)relid[2]=1 ;
512     if(relid[3]<1)relid[3]=1 ;
513     if(relid[2]>GetNPhi())relid[2]=GetNPhi() ;
514     if(relid[3]>GetNZ())relid[3]=GetNZ() ;
515     RelToAbsNumbering(relid,absId) ;
516   }
517   else{
518     Int_t indx2 = cpath.Index("/",indx) ;
519     if(indx2==-1)
520       indx2=cpath.Length() ;
521     TString cell=cpath(indx+5,indx2-indx-5) ;
522     Int_t icell=cell.Atoi() ;
523     indx = cpath.Index("PSTR") ;
524     indx2 = cpath.Index("/",indx) ;
525     TString strip=cpath(indx+5,indx2-indx-5) ;
526     Int_t iStrip = strip.Atoi() ; 
527
528     Int_t row = fGeometryEMCA->GetNStripZ() - (iStrip - 1) % (fGeometryEMCA->GetNStripZ()) ;
529     Int_t col = (Int_t) TMath::Ceil((Double_t) iStrip/(fGeometryEMCA->GetNStripZ())) -1 ;
530  
531     // Absid for 8x2-strips. Looks nice :)
532     absId = (module-1)*GetNCristalsInModule() +
533                   row * 2 + (col*fGeometryEMCA->GetNCellsXInStrip() + (icell - 1) / 2)*GetNZ() - (icell & 1 ? 1 : 0);
534  
535   }
536  
537 }
538
539 //____________________________________________________________________________
540 void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) const 
541 {
542   // Converts the relative numbering into the local PHOS-module (x, z) coordinates
543   // Note: sign of z differs from that in the previous version (Yu.Kharlov, 12 Oct 2000)
544   
545
546   if (!gGeoManager){
547     AliFatal("Geo manager not initialized\n");
548   }
549   //construct module name
550   char path[100] ;
551   if(relid[1]==0){ //this is PHOS
552
553 //   Calculations using ideal geometry (obsolete)
554 //    x = - ( GetNPhi()/2. - relid[2]    + 0.5 ) *  GetCellStep() ; // position of Xtal with respect
555 //    z = - ( GetNZ()  /2. - relid[3] + 0.5 ) *  GetCellStep() ; // of center of PHOS module  
556
557     Double_t pos[3]= {0.0,-fCryCellShift,0.}; //Position incide the crystal 
558     Double_t posC[3]={0.0,0.0,0.}; //Global position
559
560     //Shift and possibly apply misalignment corrections
561     Int_t nCellsXInStrip=fGeometryEMCA->GetNCellsXInStrip() ;
562     Int_t nCellsZInStrip=fGeometryEMCA->GetNCellsZInStrip() ;
563 //    Int_t strip=1+(relid[3]-1)/fGeometryEMCA->GetNCellsZInStrip()+((relid[2]-1)/nCellsInStrip)*fGeometryEMCA->GetNStripZ() ;
564     Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/nCellsXInStrip))*fGeometryEMCA->GetNStripZ()-
565                 (Int_t) TMath::Ceil((Double_t)relid[3]/nCellsZInStrip) ;
566     Int_t cellraw= relid[3]%nCellsZInStrip ;
567     if(cellraw==0)cellraw=nCellsZInStrip ;
568     Int_t cell= ((relid[2]-1)%nCellsXInStrip)*nCellsZInStrip + cellraw ; 
569     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1/PSTR_%d/PCEL_%d",
570             relid[0],strip,cell) ;
571     if (!gGeoManager->cd(path)){
572       AliFatal("Geo manager can not find path \n");
573     }
574     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
575     if (m) m->LocalToMaster(pos,posC);
576     else{
577       AliFatal("Geo matrixes are not loaded \n") ;
578     }
579     //Return to PHOS local system  
580     Double_t posL[3]={posC[0],posC[1],posC[2]};
581     sprintf(path,"/ALIC_1/PHOS_%d",relid[0]) ;
582     if (!gGeoManager->cd(path)){
583       AliFatal("Geo manager can not find path \n");
584     }
585     TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
586     if (mPHOS) mPHOS->MasterToLocal(posC,posL);
587     else{
588       AliFatal("Geo matrixes are not loaded \n") ;
589     }
590 //printf("RelPosInMod: posL=[%f,%f,%f]\n",posL[0],posL[1],posL[2]) ;
591 //printf("old: x=%f, z=%f \n",x,z);
592     x=posL[0] ;
593     z=posL[1];
594     return ;
595   }
596   else{//CPV
597     //first calculate position with respect to CPV plain 
598     Int_t row        = relid[2] ; //offset along x axis
599     Int_t column     = relid[3] ; //offset along z axis
600     Double_t pos[3]= {0.0,0.0,0.}; //Position incide the CPV printed circuit
601     Double_t posC[3]={0.0,0.0,0.}; //Global position
602     pos[0] = - ( GetNumberOfCPVPadsPhi()/2. - row    - 0.5 ) * GetPadSizePhi()  ; // position of pad  with respect
603     pos[2] = - ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module
604  
605     //now apply possible shifts and rotations
606     sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",relid[0]) ;
607     if (!gGeoManager->cd(path)){
608       AliFatal("Geo manager can not find path \n");
609     }
610     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
611     if (m) m->LocalToMaster(pos,posC);
612     else{
613       AliFatal("Geo matrixes are not loaded \n") ;
614     }
615     //Return to PHOS local system
616     Double_t posL[3]={0.,0.,0.,} ;
617     sprintf(path,"/ALIC_1/PHOS_%d",relid[0]) ;
618     if (!gGeoManager->cd(path)){
619       AliFatal("Geo manager can not find path \n");
620     }
621     TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
622     if (mPHOS) mPHOS->MasterToLocal(posC,posL);
623     else{
624       AliFatal("Geo matrixes are not loaded \n") ;
625     }
626     x=posL[0] ;
627     z=posL[1];
628     return ;
629  
630   }
631   
632 }
633
634 //____________________________________________________________________________
635
636 void AliPHOSGeometry::GetModuleCenter(TVector3& center, 
637                                       const char *det,
638                                       Int_t module) const
639 {
640   // Returns a position of the center of the CPV or EMC module
641   // in ideal (not misaligned) geometry
642   Float_t rDet = 0.;
643   if      (strcmp(det,"CPV") == 0) rDet  = GetIPtoCPVDistance   ();
644   else if (strcmp(det,"EMC") == 0) rDet  = GetIPtoCrystalSurface();
645   else 
646     AliFatal(Form("Wrong detector name %s",det));
647
648   Float_t angle = GetPHOSAngle(module); // (40,20,0,-20,-40) degrees
649   angle *= TMath::Pi()/180;
650   angle += 3*TMath::Pi()/2.;
651   center.SetXYZ(rDet*TMath::Cos(angle), rDet*TMath::Sin(angle), 0.);
652 }
653
654 //____________________________________________________________________________
655
656 void AliPHOSGeometry::Global2Local(TVector3& localPosition,
657                                    const TVector3& globalPosition,
658                                    Int_t module) const
659 {
660   // Transforms a global position of the rec.point to the local coordinate system
661   //Return to PHOS local system
662   Double_t posG[3]={globalPosition.X(),globalPosition.Y(),globalPosition.Z()} ;
663   Double_t posL[3]={0.,0.,0.} ;
664   char path[100] ;
665   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",module) ;
666 //  sprintf(path,"/ALIC_1/PHOS_%d",module) ;
667   if (!gGeoManager->cd(path)){
668     AliFatal("Geo manager can not find path \n");
669   }
670   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
671   if (mPHOS) mPHOS->MasterToLocal(posG,posL);
672   else{
673     AliFatal("Geo matrixes are not loaded \n") ;
674   }
675   localPosition.SetXYZ(posL[0],posL[1]+fCrystalShift,-posL[2]) ;  
676  
677 /*
678   Float_t angle = GetPHOSAngle(module); // (40,20,0,-20,-40) degrees
679   angle *= TMath::Pi()/180;
680   angle += 3*TMath::Pi()/2.;
681   localPosition = globalPosition;
682   localPosition.RotateZ(-angle);
683 */
684 }
685 //____________________________________________________________________________
686 void AliPHOSGeometry::Local2Global(Int_t mod, Float_t x, Float_t z,
687                                    TVector3& globalPosition) const 
688 {
689   char path[100] ;
690   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",mod) ;
691 //  sprintf(path,"/ALIC_1/PHOS_%d",mod) ;
692   if (!gGeoManager->cd(path)){
693     AliFatal("Geo manager can not find path \n");
694   }
695   Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!!
696   Double_t posG[3] ;
697   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
698   if (mPHOS){
699      mPHOS->LocalToMaster(posL,posG);
700   }    
701   else{
702     AliFatal("Geo matrixes are not loaded \n") ;
703   }
704   globalPosition.SetXYZ(posG[0],posG[1],posG[2]) ;
705 }
706 //____________________________________________________________________________
707 void AliPHOSGeometry::GetIncidentVector(const TVector3 &vtx, Int_t module, Float_t x,Float_t z, TVector3 &vInc) const {
708   //Calculates vector pointing from vertex to current poisition in module local frame
709   //Note that PHOS local system and ALICE global have opposite z directions
710
711   Global2Local(vInc,vtx,module) ; 
712   vInc.SetXYZ(vInc.X()+x,vInc.Y(),vInc.Z()+z) ;
713 }