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