Coding rule violations fixed.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSGeoUtils.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: AliPHOSGeometry.cxx 25590 2008-05-06 07:09:11Z prsnko $ */
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 "TParticle.h"
33 #include <TGeoManager.h>
34 #include <TGeoMatrix.h>
35
36 // --- Standard library ---
37
38 // --- AliRoot header files ---
39 #include "AliPHOSEMCAGeometry.h"
40 #include "AliPHOSCPVGeometry.h"
41 #include "AliPHOSSupportGeometry.h"
42 #include "AliPHOSGeoUtils.h"
43
44 ClassImp(AliPHOSGeoUtils)
45
46 //____________________________________________________________________________
47 AliPHOSGeoUtils::AliPHOSGeoUtils():
48   fGeometryEMCA(0x0),fGeometryCPV(0x0),fGeometrySUPP(0x0),
49   fNModules(0),fNCristalsInModule(0),fNPhi(0),fNZ(0),
50   fNumberOfCPVPadsPhi(0),fNumberOfCPVPadsZ(0),
51   fNCellsXInStrip(0),fNCellsZInStrip(0),fNStripZ(0),
52   fCrystalShift(0.),fCryCellShift(0.),fCellStep(0.),
53   fPadSizePhi(0.),fPadSizeZ(0.),fCPVBoxSizeY(0.)
54  
55 {
56     // default ctor 
57     // must be kept public for root persistency purposes, but should never be called by the outside world
58 }  
59
60 //____________________________________________________________________________
61 AliPHOSGeoUtils::AliPHOSGeoUtils(const AliPHOSGeoUtils & rhs)
62   : TNamed(rhs),
63   fGeometryEMCA(0x0),fGeometryCPV(0x0),fGeometrySUPP(0x0),
64   fNModules(0),fNCristalsInModule(0),fNPhi(0),fNZ(0),
65   fNumberOfCPVPadsPhi(0),fNumberOfCPVPadsZ(0),
66   fNCellsXInStrip(0),fNCellsZInStrip(0),fNStripZ(0),
67   fCrystalShift(0.),fCryCellShift(0.),fCellStep(0.),
68   fPadSizePhi(0.),fPadSizeZ(0.),fCPVBoxSizeY(0.)
69 {
70   Fatal("cpy ctor", "not implemented") ; 
71 }
72
73 //____________________________________________________________________________
74 AliPHOSGeoUtils::AliPHOSGeoUtils(const Text_t* name, const Text_t* title) 
75     : TNamed(name, title),
76   fGeometryEMCA(0x0),fGeometryCPV(0x0),fGeometrySUPP(0x0),
77   fNModules(0),fNCristalsInModule(0),fNPhi(0),fNZ(0),
78   fNumberOfCPVPadsPhi(0),fNumberOfCPVPadsZ(0),
79   fNCellsXInStrip(0),fNCellsZInStrip(0),fNStripZ(0),
80   fCrystalShift(0.),fCryCellShift(0.),fCellStep(0.),
81   fPadSizePhi(0.),fPadSizeZ(0.),fCPVBoxSizeY(0.)
82
83   // ctor only for normal usage 
84
85   fGeometryEMCA = new AliPHOSEMCAGeometry() ;
86   fGeometryCPV  = new AliPHOSCPVGeometry() ;
87   fGeometrySUPP = new AliPHOSSupportGeometry() ;
88
89   fNModules     = 5;
90   fNPhi  = fGeometryEMCA->GetNPhi() ;
91   fNZ    = fGeometryEMCA->GetNZ() ;
92   fNCristalsInModule = fNPhi*fNZ ;
93   fNCellsXInStrip= fGeometryEMCA->GetNCellsXInStrip() ;
94   fNCellsZInStrip= fGeometryEMCA->GetNCellsZInStrip() ;
95   fNStripZ = fGeometryEMCA->GetNStripZ() ;
96   fXtlArrSize[0]=fGeometryEMCA->GetInnerThermoHalfSize()[0] ; //Wery close to the zise of the Xtl set
97   fXtlArrSize[1]=fGeometryEMCA->GetInnerThermoHalfSize()[1] ; //Wery close to the zise of the Xtl set
98   fXtlArrSize[2]=fGeometryEMCA->GetInnerThermoHalfSize()[2] ; //Wery close to the zise of the Xtl set
99
100   //calculate offset to crystal surface
101   Float_t * inthermo = fGeometryEMCA->GetInnerThermoHalfSize() ;
102   Float_t * strip = fGeometryEMCA->GetStripHalfSize() ;
103   Float_t* splate = fGeometryEMCA->GetSupportPlateHalfSize();
104   Float_t * crystal = fGeometryEMCA->GetCrystalHalfSize() ;
105   Float_t * pin = fGeometryEMCA->GetAPDHalfSize() ;
106   Float_t * preamp = fGeometryEMCA->GetPreampHalfSize() ;
107   fCrystalShift=-inthermo[1]+strip[1]+splate[1]+crystal[1]-fGeometryEMCA->GetAirGapLed()/2.+pin[1]+preamp[1] ;
108   fCryCellShift=crystal[1]-(fGeometryEMCA->GetAirGapLed()-2*pin[1]-2*preamp[1])/2;
109   fCellStep = 2.*fGeometryEMCA->GetAirCellHalfSize()[0] ;
110
111
112   fNumberOfCPVPadsPhi = fGeometryCPV->GetNumberOfCPVPadsPhi() ;
113   fNumberOfCPVPadsZ   = fGeometryCPV->GetNumberOfCPVPadsZ() ;
114   fPadSizePhi = fGeometryCPV->GetCPVPadSizePhi() ;
115   fPadSizeZ   = fGeometryCPV->GetCPVPadSizeZ() ; 
116   fCPVBoxSizeY= fGeometryCPV->GetCPVBoxSize(1) ;
117 }
118
119 //____________________________________________________________________________
120 AliPHOSGeoUtils & AliPHOSGeoUtils::operator = (const AliPHOSGeoUtils  & /*rvalue*/) { 
121
122   Fatal("assignment operator", "not implemented") ; 
123     return *this ;
124 }
125
126 //____________________________________________________________________________
127 AliPHOSGeoUtils::~AliPHOSGeoUtils(void)
128 {
129   // dtor
130   if(fGeometryEMCA){
131     delete fGeometryEMCA; fGeometryEMCA = 0 ;
132   }
133   if(fGeometryCPV){
134     delete fGeometryCPV; fGeometryCPV=0 ;
135   }
136   if(fGeometrySUPP){
137     delete fGeometrySUPP ; fGeometrySUPP=0 ;
138   }
139
140 }
141 //____________________________________________________________________________
142 Bool_t AliPHOSGeoUtils::AbsToRelNumbering(Int_t absId, Int_t * relid) const
143 {
144   // Converts the absolute numbering into the following array
145   //  relid[0] = PHOS Module number 1:fNModules 
146   //  relid[1] = 0 if PbW04
147   //           = -1 if CPV
148   //  relid[2] = Row number inside a PHOS module
149   //  relid[3] = Column number inside a PHOS module
150
151   Float_t id = absId ;
152
153   Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / fNCristalsInModule ) ; 
154   
155   if ( phosmodulenumber >  fNModules ) { // it is a CPV pad
156     
157     id -=  fNPhi * fNZ *  fNModules ; 
158     Float_t nCPV  = fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ ;
159     relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ;
160     relid[1] = -1 ;
161     id -= ( relid[0] - 1 ) * nCPV ; 
162     relid[2] = (Int_t) TMath::Ceil( id / fNumberOfCPVPadsZ ) ;
163     relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * fNumberOfCPVPadsZ ) ; 
164   } 
165   else { // it is a PW04 crystal
166
167     relid[0] = phosmodulenumber ;
168     relid[1] = 0 ;
169     id -= ( phosmodulenumber - 1 ) *  fNPhi * fNZ ; 
170     relid[2] = (Int_t)TMath::Ceil( id / fNZ )  ;
171     relid[3] = (Int_t)( id - ( relid[2] - 1 ) * fNZ ) ; 
172   } 
173   return kTRUE ; 
174 }
175 //____________________________________________________________________________
176 Bool_t AliPHOSGeoUtils::RelToAbsNumbering(const Int_t * relid, Int_t &  absId) const
177 {
178   // Converts the relative numbering into the absolute numbering
179   // EMCA crystals:
180   //  absId = from 1 to fNModules * fNPhi * fNZ
181   // CPV pad:
182   //  absId = from N(total PHOS crystals) + 1
183   //          to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ
184
185   if ( relid[1] ==  0 ) {                            // it is a Phos crystal
186     absId =
187       ( relid[0] - 1 ) * fNPhi * fNZ         // the offset of PHOS modules
188       + ( relid[2] - 1 ) * fNZ                   // the offset along phi
189       +   relid[3] ;                                 // the offset along z
190   }
191   else { // it is a CPV pad
192     absId =    fNPhi * fNZ *  fNModules         // the offset to separate EMCA crystals from CPV pads
193       + ( relid[0] - 1 ) * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ   // the pads offset of PHOS modules 
194       + ( relid[2] - 1 ) * fNumberOfCPVPadsZ                            // the pads offset of a CPV row
195       +   relid[3] ;                                                         // the column number
196   }
197   
198   return kTRUE ; 
199 }
200
201 //____________________________________________________________________________
202 void AliPHOSGeoUtils::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) const 
203 {
204   // Converts the relative numbering into the local PHOS-module (x, z) coordinates
205
206   if (!gGeoManager){
207     printf("Geo manager not initialized\n");
208     abort() ;
209   }
210   //construct module name
211   char path[100] ;
212   if(relid[1]==0){ //this is PHOS
213
214     Double_t pos[3]= {0.0,-fCryCellShift,0.}; //Position incide the crystal 
215     Double_t posC[3]={0.0,0.0,0.}; //Global position
216
217     //Shift and possibly apply misalignment corrections
218     Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/fNCellsXInStrip))*fNStripZ-
219                 (Int_t) TMath::Ceil((Double_t)relid[3]/fNCellsZInStrip) ;
220     Int_t cellraw= relid[3]%fNCellsZInStrip ;
221     if(cellraw==0)cellraw=fNCellsZInStrip ;
222     Int_t cell= ((relid[2]-1)%fNCellsXInStrip)*fNCellsZInStrip + cellraw ; 
223     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1/PSTR_%d/PCEL_%d",
224             relid[0],strip,cell) ;
225     if (!gGeoManager->cd(path)){
226       printf("Geo manager can not find path \n");
227       abort() ;
228     }
229     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
230     if (m) m->LocalToMaster(pos,posC);
231     else{
232       printf("Geo matrixes are not loaded \n") ;
233       abort() ;
234     }
235     //Return to PHOS local system  
236     Double_t posL[3]={posC[0],posC[1],posC[2]};
237     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",relid[0]) ;
238     //    sprintf(path,"/ALIC_1/PHOS_%d",relid[0]) ;
239     if (!gGeoManager->cd(path)){
240       printf("Geo manager can not find path \n");
241       abort();
242     }
243     TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
244     if (mPHOS) mPHOS->MasterToLocal(posC,posL);
245     else{
246       printf("Geo matrixes are not loaded \n") ;
247       abort() ;
248     }
249     x=posL[0] ;
250     z=-posL[2];
251     return ;
252   }
253   else{//CPV
254     //first calculate position with respect to CPV plain 
255     Int_t row        = relid[2] ; //offset along x axis
256     Int_t column     = relid[3] ; //offset along z axis
257     Double_t pos[3]= {0.0,0.0,0.}; //Position incide the CPV printed circuit
258     Double_t posC[3]={0.0,0.0,0.}; //Global position
259     pos[0] = - ( fNumberOfCPVPadsPhi/2. - row    - 0.5 ) * fPadSizePhi  ; // position of pad  with respect
260     pos[2] = - ( fNumberOfCPVPadsZ  /2. - column - 0.5 ) * fPadSizeZ  ; // of center of PHOS module
261
262     //now apply possible shifts and rotations
263     sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",relid[0]) ;
264     if (!gGeoManager->cd(path)){
265       printf("Geo manager can not find path \n");
266       abort() ;
267     }
268     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
269     if (m) m->LocalToMaster(pos,posC);
270     else{
271       printf("Geo matrixes are not loaded \n") ;
272       abort() ;
273     }
274     //Return to PHOS local system
275     Double_t posL[3]={0.,0.,0.,} ;
276     sprintf(path,"/ALIC_1/PHOS_%d",relid[0]) ;
277     if (!gGeoManager->cd(path)){
278       printf("Geo manager can not find path \n");
279       abort() ;
280     }
281     TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
282     if (mPHOS) mPHOS->MasterToLocal(posC,posL);
283     else{
284       printf("Geo matrixes are not loaded \n") ;
285       abort() ;
286     }
287     x=posL[0] ;
288     z=posL[1];
289     return ;
290  
291   }
292   
293 }
294 //____________________________________________________________________________
295 void AliPHOSGeoUtils::RelPosToAbsId(Int_t module, Double_t x, Double_t z, Int_t & absId) const
296 {
297   // converts local PHOS-module (x, z) coordinates to absId 
298
299   //find Global position
300   if (!gGeoManager){
301     printf("Geo manager not initialized\n");
302     abort() ;
303   }
304   Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!!
305   Double_t posG[3] ;
306   char path[100] ;
307   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",module) ;
308   if (!gGeoManager->cd(path)){
309     printf("Geo manager can not find path \n");
310     abort() ;
311   }
312   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
313   if (mPHOS){
314      mPHOS->LocalToMaster(posL,posG);
315   }
316   else{
317     printf("Geo matrixes are not loaded \n") ;
318     abort() ;
319   }
320
321   Int_t relid[4] ;
322   gGeoManager->FindNode(posG[0],posG[1],posG[2]) ;
323   //Check that path contains PSTR and extract strip number
324   TString cpath(gGeoManager->GetPath()) ;
325   Int_t indx = cpath.Index("PCEL") ;
326   if(indx==-1){ //for the few events when particle hits between srips use ideal geometry
327     relid[0] = module ;
328     relid[1] = 0 ;
329     relid[2] = static_cast<Int_t>(TMath::Ceil( x/ fCellStep + fNPhi / 2.) );
330     relid[3] = static_cast<Int_t>(TMath::Ceil(-z/ fCellStep + fNZ   / 2.) ) ;
331     if(relid[2]<1)relid[2]=1 ;
332     if(relid[3]<1)relid[3]=1 ;
333     if(relid[2]>fNPhi)relid[2]=fNPhi ;
334     if(relid[3]>fNZ)relid[3]=fNZ ;
335     RelToAbsNumbering(relid,absId) ;
336   }
337   else{
338     Int_t indx2 = cpath.Index("/",indx) ;
339     if(indx2==-1)
340       indx2=cpath.Length() ;
341     TString cell=cpath(indx+5,indx2-indx-5) ;
342     Int_t icell=cell.Atoi() ;
343     indx = cpath.Index("PSTR") ;
344     indx2 = cpath.Index("/",indx) ;
345     TString strip=cpath(indx+5,indx2-indx-5) ;
346     Int_t iStrip = strip.Atoi() ; 
347
348     Int_t row = fNStripZ - (iStrip - 1) % (fNStripZ) ;
349     Int_t col = (Int_t) TMath::Ceil((Double_t) iStrip/(fNStripZ)) -1 ;
350  
351     // Absid for 8x2-strips. Looks nice :)
352     absId = (module-1)*fNCristalsInModule +
353                   row * 2 + (col*fNCellsXInStrip + (icell - 1) / 2)*fNZ - (icell & 1 ? 1 : 0);
354  
355   }
356  
357 }
358
359 //____________________________________________________________________________
360 void AliPHOSGeoUtils::RelPosToRelId(Int_t module, Double_t x, Double_t z, Int_t * relId) const
361 {
362   //Evaluates RelId of the crystall with given coordinates
363
364   Int_t absId ;
365   RelPosToAbsId(module, x,z,absId) ;
366   AbsToRelNumbering(absId,relId) ;
367 }
368
369 //____________________________________________________________________________
370 void AliPHOSGeoUtils::RelPosInAlice(Int_t id, TVector3 & pos ) const
371 {
372   // Converts the absolute numbering into the global ALICE coordinate system
373   
374   if (!gGeoManager){
375     printf("Geo manager not initialized\n");
376     abort();
377   }
378     
379   Int_t relid[4] ;
380     
381   AbsToRelNumbering(id , relid) ;
382     
383   //construct module name
384   char path[100] ;
385   if(relid[1]==0){ //this is EMC
386  
387     Double_t ps[3]= {0.0,-fCryCellShift,0.}; //Position incide the crystal 
388     Double_t psC[3]={0.0,0.0,0.}; //Global position
389  
390     //Shift and possibly apply misalignment corrections
391     Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/fNCellsXInStrip))*fNStripZ-
392                 (Int_t) TMath::Ceil((Double_t)relid[3]/fNCellsZInStrip) ;
393     Int_t cellraw= relid[3]%fNCellsZInStrip ;
394     if(cellraw==0)cellraw=fNCellsZInStrip ;
395     Int_t cell= ((relid[2]-1)%fNCellsXInStrip)*fNCellsZInStrip + cellraw ;
396     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1/PSTR_%d/PCEL_%d",
397             relid[0],strip,cell) ;
398     if (!gGeoManager->cd(path)){
399       printf("Geo manager can not find path \n");
400       abort() ;
401     }
402     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
403     if (m) m->LocalToMaster(ps,psC);
404     else{
405       printf("Geo matrixes are not loaded \n") ;
406       abort() ;
407     }
408     pos.SetXYZ(psC[0],psC[1],psC[2]) ; 
409   }
410   else{
411     //first calculate position with respect to CPV plain
412     Int_t row        = relid[2] ; //offset along x axis
413     Int_t column     = relid[3] ; //offset along z axis
414     Double_t ps[3]= {0.0,fCPVBoxSizeY/2.,0.}; //Position on top of CPV
415     Double_t psC[3]={0.0,0.0,0.}; //Global position
416     pos[0] = - ( fNumberOfCPVPadsPhi/2. - row    - 0.5 ) * fPadSizePhi  ; // position of pad  with respect
417     pos[2] = - ( fNumberOfCPVPadsZ  /2. - column - 0.5 ) * fPadSizeZ  ; // of center of PHOS module
418  
419     //now apply possible shifts and rotations
420     sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",relid[0]) ;
421     if (!gGeoManager->cd(path)){
422       printf("Geo manager can not find path \n");
423       abort();
424     }
425     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
426     if (m) m->LocalToMaster(ps,psC);
427     else{
428       printf("Geo matrixes are not loaded \n") ;
429       abort() ;
430     }
431     pos.SetXYZ(psC[0],psC[1],-psC[2]) ; 
432   }
433
434
435 //____________________________________________________________________________
436 void AliPHOSGeoUtils::Local2Global(Int_t mod, Float_t x, Float_t z,
437                                    TVector3& globalPosition) const 
438 {
439   char path[100] ;
440   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",mod) ;
441   if (!gGeoManager->cd(path)){
442     printf("Geo manager can not find path \n");
443     abort() ;
444   }
445   Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!!
446   Double_t posG[3] ;
447   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
448   if (mPHOS){
449      mPHOS->LocalToMaster(posL,posG);
450   }    
451   else{
452     printf("Geo matrixes are not loaded \n") ;
453     abort() ;
454   }
455   globalPosition.SetXYZ(posG[0],posG[1],posG[2]) ;
456 }
457 //____________________________________________________________________________
458 void AliPHOSGeoUtils::Global2Local(TVector3& localPosition,
459                                    const TVector3& globalPosition,
460                                    Int_t module) const
461 {
462   // Transforms a global position to the local coordinate system
463   // of the module 
464   //Return to PHOS local system
465   Double_t posG[3]={globalPosition.X(),globalPosition.Y(),globalPosition.Z()} ;
466   Double_t posL[3]={0.,0.,0.} ;
467   char path[100] ;
468   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",module) ;
469   if (!gGeoManager->cd(path)){
470     printf("Geo manager can not find path \n");
471     abort() ;
472   }
473   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
474   if (mPHOS) mPHOS->MasterToLocal(posG,posL);
475   else{
476     printf("Geo matrixes are not loaded \n") ;
477     abort() ;
478   }
479   localPosition.SetXYZ(posL[0],posL[1]+fCrystalShift,-posL[2]) ;  
480  
481 }
482 //____________________________________________________________________________
483 Bool_t AliPHOSGeoUtils::GlobalPos2RelId(TVector3 & global, Int_t * relId){
484   //Converts position in global ALICE coordinates to relId 
485   //returns false if x,z coordinates are beyond PHOS
486   //distande to PHOS surface is NOT calculated 
487   TVector3 loc ;
488   for(Int_t mod=1; mod<fNModules; mod++){
489     Global2Local(loc,global,mod) ;
490     //If in Acceptance
491     if((TMath::Abs(loc.Z())<fXtlArrSize[2]) && (TMath::Abs(loc.X())<fXtlArrSize[0])){
492        RelPosToRelId(mod,loc.X(),loc.Z(),relId);
493        return kTRUE ;
494     }
495   }
496   return kFALSE ; 
497
498 }
499 //____________________________________________________________________________
500 Bool_t AliPHOSGeoUtils::ImpactOnEmc(const TParticle * particle,
501        Int_t & moduleNumber, Double_t & z, Double_t & x) const
502 {
503   // Tells if a particle enters PHOS and evaluates hit position
504   Double_t vtx[3]={particle->Vx(),particle->Vy(),particle->Vz()} ;
505   return ImpactOnEmc(vtx,particle->Theta(),particle->Phi(),moduleNumber,z,x);
506 }
507  
508 //____________________________________________________________________________
509 Bool_t AliPHOSGeoUtils::ImpactOnEmc(const Double_t * vtx, Double_t theta, Double_t phi, 
510                                   Int_t & moduleNumber, Double_t & z, Double_t & x) const
511 {
512   // calculates the impact coordinates on PHOS of a neutral particle  
513   // emitted in the vertex vtx[3] with direction vec(p) in the ALICE global coordinate system
514   TVector3 p(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) ;
515   return ImpactOnEmc(vtx,p,moduleNumber,z,x) ;
516
517 }
518 //____________________________________________________________________________
519 Bool_t AliPHOSGeoUtils::ImpactOnEmc(const Double_t * vtx, const TVector3 &p,
520                                   Int_t & moduleNumber, Double_t & z, Double_t & x) const
521 {
522   // calculates the impact coordinates on PHOS of a neutral particle  
523   // emitted in the vertex vtx[3] with direction theta and phi in the ALICE global coordinate system
524   TVector3 v(vtx[0],vtx[1],vtx[2]) ;
525
526   if (!gGeoManager){
527     printf("Geo manager not initialized\n");
528     abort() ;
529     return kFALSE ;
530   }
531  
532   for(Int_t imod=1; imod<=fNModules ; imod++){
533     //create vector from (0,0,0) to center of crystal surface of imod module
534     Double_t tmp[3]={0.,-fCrystalShift,0.} ;
535  
536     char path[100] ;
537     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",imod) ;
538     if (!gGeoManager->cd(path)){
539       printf("Geo manager can not find path \n");
540       abort() ;
541       return kFALSE ;
542     }
543     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
544     Double_t posG[3]={0.,0.,0.} ;
545     if (m) m->LocalToMaster(tmp,posG);
546     TVector3 n(posG[0],posG[1],posG[2]) ; 
547     Double_t direction=n.Dot(p) ;
548     if(direction<=0.)
549       continue ; //momentum directed FROM module
550     Double_t fr = (n.Mag2()-n.Dot(v))/direction ;  
551     //Calculate direction in module plain
552     n-=v+fr*p ;
553     n*=-1. ;
554     if(TMath::Abs(TMath::Abs(n.Z())<fXtlArrSize[2]) && n.Pt()<fXtlArrSize[0]){
555       moduleNumber = imod ;
556       z=n.Z() ;
557       x=TMath::Sign(n.Pt(),n.X()) ;
558       //no need to return to local system since we calcilated distance from module center
559       //and tilts can not be significant.
560       return kTRUE ;
561     }
562   }
563   //Not in acceptance
564   x=0; z=0 ;
565   moduleNumber=0 ;
566   return kFALSE ;
567
568 }
569 //____________________________________________________________________________
570 void AliPHOSGeoUtils::GetIncidentVector(const TVector3 &vtx, Int_t module, Float_t x,Float_t z, TVector3 &vInc) const {
571   //Calculates vector pointing from vertex to current poisition in module local frame
572   //Note that PHOS local system and ALICE global have opposite z directions
573
574   Global2Local(vInc,vtx,module) ; 
575   vInc.SetXYZ(vInc.X()+x,vInc.Y(),vInc.Z()+z) ;
576 }
577