7d70510ee8e375cdfc39446bc2fb6fd496ebd869
[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 
20 // PHOS consists of the electromagnetic calorimeter (EMCA)
21 // and a charged particle veto (CPV)
22 // The EMCA/CPV modules are parametrized so that any configuration
23 // can be easily implemented 
24 // The title is used to identify the version of CPV used.
25 //                  
26 // -- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (RRC "KI" & SUBATECH)
27
28 // --- ROOT system ---
29
30 #include "TClonesArray.h"
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 "AliLog.h"
40 #include "AliPHOSEMCAGeometry.h"
41 #include "AliPHOSCPVGeometry.h"
42 #include "AliPHOSSupportGeometry.h"
43 #include "AliPHOSGeoUtils.h"
44
45 ClassImp(AliPHOSGeoUtils)
46
47 //____________________________________________________________________________
48 AliPHOSGeoUtils::AliPHOSGeoUtils():
49   fGeometryEMCA(0x0),fGeometryCPV(0x0),fGeometrySUPP(0x0),
50   fNModules(0),fNCristalsInModule(0),fNPhi(0),fNZ(0),
51   fNumberOfCPVPadsPhi(0),fNumberOfCPVPadsZ(0),
52   fNCellsXInStrip(0),fNCellsZInStrip(0),fNStripZ(0),
53   fCrystalShift(0.),fCryCellShift(0.),fCryStripShift(0.),fCellStep(0.),
54   fPadSizePhi(0.),fPadSizeZ(0.),fCPVBoxSizeY(0.),fMisalArray(0x0)
55  
56 {
57     // default ctor 
58     // must be kept public for root persistency purposes, but should never be called by the outside world
59 }  
60
61 //____________________________________________________________________________
62 AliPHOSGeoUtils::AliPHOSGeoUtils(const AliPHOSGeoUtils & rhs)
63   : TNamed(rhs),
64   fGeometryEMCA(0x0),fGeometryCPV(0x0),fGeometrySUPP(0x0),
65   fNModules(0),fNCristalsInModule(0),fNPhi(0),fNZ(0),
66   fNumberOfCPVPadsPhi(0),fNumberOfCPVPadsZ(0),
67   fNCellsXInStrip(0),fNCellsZInStrip(0),fNStripZ(0),
68   fCrystalShift(0.),fCryCellShift(0.),fCryStripShift(0.),fCellStep(0.),
69   fPadSizePhi(0.),fPadSizeZ(0.),fCPVBoxSizeY(0.),fMisalArray(0x0)
70 {
71   Fatal("cpy ctor", "not implemented") ; 
72 }
73
74 //____________________________________________________________________________
75 AliPHOSGeoUtils::AliPHOSGeoUtils(const Text_t* name, const Text_t* title) 
76     : TNamed(name, title),
77   fGeometryEMCA(0x0),fGeometryCPV(0x0),fGeometrySUPP(0x0),
78   fNModules(0),fNCristalsInModule(0),fNPhi(0),fNZ(0),
79   fNumberOfCPVPadsPhi(0),fNumberOfCPVPadsZ(0),
80   fNCellsXInStrip(0),fNCellsZInStrip(0),fNStripZ(0),
81   fCrystalShift(0.),fCryCellShift(0.),fCryStripShift(0.),fCellStep(0.),
82   fPadSizePhi(0.),fPadSizeZ(0.),fCPVBoxSizeY(0.),fMisalArray(0x0)
83
84   // ctor only for normal usage 
85
86   fGeometryEMCA = new AliPHOSEMCAGeometry() ;
87   fGeometryCPV  = new AliPHOSCPVGeometry() ;
88   fGeometrySUPP = new AliPHOSSupportGeometry() ;
89
90   fNModules     = 5;
91   fNPhi  = fGeometryEMCA->GetNPhi() ;
92   fNZ    = fGeometryEMCA->GetNZ() ;
93   fNCristalsInModule = fNPhi*fNZ ;
94   fNCellsXInStrip= fGeometryEMCA->GetNCellsXInStrip() ;
95   fNCellsZInStrip= fGeometryEMCA->GetNCellsZInStrip() ;
96   fNStripZ = fGeometryEMCA->GetNStripZ() ;
97   fXtlArrSize[0]=fGeometryEMCA->GetInnerThermoHalfSize()[0] ; //Wery close to the zise of the Xtl set
98   fXtlArrSize[1]=fGeometryEMCA->GetInnerThermoHalfSize()[1] ; //Wery close to the zise of the Xtl set
99   fXtlArrSize[2]=fGeometryEMCA->GetInnerThermoHalfSize()[2] ; //Wery close to the zise of the Xtl set
100
101   //calculate offset to crystal surface
102   const Float_t * inthermo = fGeometryEMCA->GetInnerThermoHalfSize() ;
103   const Float_t * strip    = fGeometryEMCA->GetStripHalfSize() ;
104   const Float_t * splate   = fGeometryEMCA->GetSupportPlateHalfSize();
105   const Float_t * crystal  = fGeometryEMCA->GetCrystalHalfSize() ;
106   const Float_t * pin      = fGeometryEMCA->GetAPDHalfSize() ;
107   const Float_t * preamp   = fGeometryEMCA->GetPreampHalfSize() ;
108   fCrystalShift=-inthermo[1]+strip[1]+splate[1]+crystal[1]-fGeometryEMCA->GetAirGapLed()/2.+pin[1]+preamp[1] ;
109   fCryCellShift=crystal[1]-(fGeometryEMCA->GetAirGapLed()-2*pin[1]-2*preamp[1])/2;
110   fCryStripShift=fCryCellShift+splate[1] ;
111   fCellStep = 2.*fGeometryEMCA->GetAirCellHalfSize()[0] ;
112
113   fNumberOfCPVPadsPhi = fGeometryCPV->GetNumberOfCPVPadsPhi() ;
114   fNumberOfCPVPadsZ   = fGeometryCPV->GetNumberOfCPVPadsZ() ;
115   fPadSizePhi = fGeometryCPV->GetCPVPadSizePhi() ;
116   fPadSizeZ   = fGeometryCPV->GetCPVPadSizeZ() ; 
117   fCPVBoxSizeY= fGeometryCPV->GetCPVBoxSize(1) ;
118
119   for(Int_t mod=0; mod<5; mod++){
120     fEMCMatrix[mod]=0 ;
121     for(Int_t istrip=0; istrip<224; istrip++)
122       fStripMatrix[mod][istrip]=0 ;
123     fCPVMatrix[mod]=0;
124     fPHOSMatrix[mod]=0 ;
125   }
126  
127 }
128
129 //____________________________________________________________________________
130 AliPHOSGeoUtils & AliPHOSGeoUtils::operator = (const AliPHOSGeoUtils  & /*rvalue*/) { 
131
132   Fatal("assignment operator", "not implemented") ; 
133     return *this ;
134 }
135
136 //____________________________________________________________________________
137 AliPHOSGeoUtils::~AliPHOSGeoUtils(void)
138 {
139   // dtor
140   if(fGeometryEMCA){
141     delete fGeometryEMCA; fGeometryEMCA = 0 ;
142   }
143   if(fGeometryCPV){
144     delete fGeometryCPV; fGeometryCPV=0 ;
145   }
146   if(fGeometrySUPP){
147     delete fGeometrySUPP ; fGeometrySUPP=0 ;
148   }
149   if(fMisalArray){
150     delete fMisalArray; fMisalArray=0 ;
151   }
152
153   delete [] fStripMatrix;
154   delete [] fEMCMatrix;
155   delete [] fCPVMatrix;
156 }
157 //____________________________________________________________________________
158 Bool_t AliPHOSGeoUtils::AbsToRelNumbering(Int_t absId, Int_t * relid) const
159 {
160   // Converts the absolute numbering into the following array
161   //  relid[0] = PHOS Module number 1:fNModules 
162   //  relid[1] = 0 if PbW04
163   //           = -1 if CPV
164   //  relid[2] = Row number inside a PHOS module
165   //  relid[3] = Column number inside a PHOS module
166
167   Float_t id = absId ;
168
169   Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / fNCristalsInModule ) ; 
170   
171   if ( phosmodulenumber >  fNModules ) { // it is a CPV pad
172     
173     id -=  fNPhi * fNZ *  fNModules ; 
174     Float_t nCPV  = fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ ;
175     relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ;
176     relid[1] = -1 ;
177     id -= ( relid[0] - 1 ) * nCPV ; 
178     relid[2] = (Int_t) TMath::Ceil( id / fNumberOfCPVPadsZ ) ;
179     relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * fNumberOfCPVPadsZ ) ; 
180   } 
181   else { // it is a PW04 crystal
182
183     relid[0] = phosmodulenumber ;
184     relid[1] = 0 ;
185     id -= ( phosmodulenumber - 1 ) *  fNPhi * fNZ ; 
186     relid[2] = (Int_t)TMath::Ceil( id / fNZ )  ;
187     relid[3] = (Int_t)( id - ( relid[2] - 1 ) * fNZ ) ; 
188   } 
189   return kTRUE ; 
190 }
191 //____________________________________________________________________________
192 Bool_t AliPHOSGeoUtils::RelToAbsNumbering(const Int_t * relid, Int_t &  absId) const
193 {
194   // Converts the relative numbering into the absolute numbering
195   // EMCA crystals:
196   //  absId = from 1 to fNModules * fNPhi * fNZ
197   // CPV pad:
198   //  absId = from N(total PHOS crystals) + 1
199   //          to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ
200
201   if ( relid[1] ==  0 ) {                            // it is a Phos crystal
202     absId =
203       ( relid[0] - 1 ) * fNPhi * fNZ         // the offset of PHOS modules
204       + ( relid[2] - 1 ) * fNZ                   // the offset along phi
205       +   relid[3] ;                                 // the offset along z
206   }
207   else { // it is a CPV pad
208     absId =    fNPhi * fNZ *  fNModules         // the offset to separate EMCA crystals from CPV pads
209       + ( relid[0] - 1 ) * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ   // the pads offset of PHOS modules 
210       + ( relid[2] - 1 ) * fNumberOfCPVPadsZ                            // the pads offset of a CPV row
211       +   relid[3] ;                                                         // the column number
212   }
213   
214   return kTRUE ; 
215 }
216
217 //____________________________________________________________________________
218 void AliPHOSGeoUtils::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) const 
219 {
220   // Converts the relative numbering into the local PHOS-module (x, z) coordinates
221
222   if(relid[1]==0){ //this is PHOS
223
224     Double_t pos[3]= {0.0,-fCryCellShift,0.}; //Position incide the crystal 
225     Double_t posC[3]={0.0,0.0,0.}; //Global position
226
227     //Shift and possibly apply misalignment corrections
228     Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/fNCellsXInStrip))*fNStripZ-
229                 (Int_t) TMath::Ceil((Double_t)relid[3]/fNCellsZInStrip) ;
230     pos[0]=((relid[2]-1)%fNCellsXInStrip-fNCellsXInStrip/2+0.5)*fCellStep ;
231     pos[2]=(-(relid[3]-1)%fNCellsZInStrip+fNCellsZInStrip/2-0.5)*fCellStep ;
232
233     Int_t mod = relid[0] ;
234     const TGeoHMatrix * m2 = GetMatrixForStrip(mod, strip) ;
235     m2->LocalToMaster(pos,posC);
236
237     //Return to PHOS local system  
238     Double_t posL2[3]={posC[0],posC[1],posC[2]};
239     const TGeoHMatrix *mPHOS2 = GetMatrixForModule(mod) ;
240     mPHOS2->MasterToLocal(posC,posL2);
241     x=posL2[0] ;
242     z=-posL2[2];
243     return ;
244   }
245   else{//CPV
246     //first calculate position with respect to CPV plain 
247     Int_t row        = relid[2] ; //offset along x axis
248     Int_t column     = relid[3] ; //offset along z axis
249     Double_t pos[3]= {0.0,0.0,0.}; //Position incide the CPV printed circuit
250     Double_t posC[3]={0.0,0.0,0.}; //Global position
251     pos[0] = - ( fNumberOfCPVPadsPhi/2. - row    - 0.5 ) * fPadSizePhi  ; // position of pad  with respect
252     pos[2] = - ( fNumberOfCPVPadsZ  /2. - column - 0.5 ) * fPadSizeZ  ; // of center of PHOS module
253
254     //now apply possible shifts and rotations
255     const TGeoHMatrix *m = GetMatrixForCPV(relid[0]) ;
256     m->LocalToMaster(pos,posC);
257     //Return to PHOS local system
258     Double_t posL[3]={0.,0.,0.,} ;
259     const TGeoHMatrix *mPHOS = GetMatrixForPHOS(relid[0]) ;
260     mPHOS->MasterToLocal(posC,posL);
261     x=posL[0] ;
262     z=posL[1];
263     return ;
264  
265   }
266   
267 }
268 //____________________________________________________________________________
269 void AliPHOSGeoUtils::RelPosToAbsId(Int_t module, Double_t x, Double_t z, Int_t & absId) const
270 {
271   // converts local PHOS-module (x, z) coordinates to absId 
272
273   //Calculate AbsId using ideal geometry. Should be sufficient for primary particles calculation
274   //(the only place where this method used currently)
275   Int_t relid[4]={module,0,1,1} ;
276   relid[2] = static_cast<Int_t>(TMath::Ceil( x/ fCellStep + fNPhi / 2.) );
277   relid[3] = fNZ+1-static_cast<Int_t>(TMath::Ceil(-z/ fCellStep + fNZ   / 2.) ) ;
278   if(relid[2]<1)relid[2]=1 ;
279   if(relid[3]<1)relid[3]=1 ;
280   if(relid[2]>fNPhi)relid[2]=fNPhi ;
281   if(relid[3]>fNZ)relid[3]=fNZ ;
282   RelToAbsNumbering(relid,absId) ;
283
284 /*
285   //find Global position
286   if (!gGeoManager){
287     printf("Geo manager not initialized\n");
288     abort() ;
289   }
290   Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!!
291   Double_t posG[3] ;
292   char path[100] ;
293   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",module) ;
294   if (!gGeoManager->cd(path)){
295     printf("Geo manager can not find path \n");
296     abort() ;
297   }
298   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
299   if (mPHOS){
300      mPHOS->LocalToMaster(posL,posG);
301   }
302   else{
303     printf("Geo matrixes are not loaded \n") ;
304     abort() ;
305   }
306
307   Int_t relid[4] ;
308   gGeoManager->FindNode(posG[0],posG[1],posG[2]) ;
309   //Check that path contains PSTR and extract strip number
310   TString cpath(gGeoManager->GetPath()) ;
311   Int_t indx = cpath.Index("PCEL") ;
312   if(indx==-1){ //for the few events when particle hits between srips use ideal geometry
313     relid[0] = module ;
314     relid[1] = 0 ;
315     relid[2] = static_cast<Int_t>(TMath::Ceil( x/ fCellStep + fNPhi / 2.) );
316     relid[3] = static_cast<Int_t>(TMath::Ceil(-z/ fCellStep + fNZ   / 2.) ) ;
317     if(relid[2]<1)relid[2]=1 ;
318     if(relid[3]<1)relid[3]=1 ;
319     if(relid[2]>fNPhi)relid[2]=fNPhi ;
320     if(relid[3]>fNZ)relid[3]=fNZ ;
321     RelToAbsNumbering(relid,absId) ;
322   }
323   else{
324     Int_t indx2 = cpath.Index("/",indx) ;
325     if(indx2==-1)
326       indx2=cpath.Length() ;
327     TString cell=cpath(indx+5,indx2-indx-5) ;
328     Int_t icell=cell.Atoi() ;
329     indx = cpath.Index("PSTR") ;
330     indx2 = cpath.Index("/",indx) ;
331     TString strip=cpath(indx+5,indx2-indx-5) ;
332     Int_t iStrip = strip.Atoi() ; 
333
334     Int_t row = fNStripZ - (iStrip - 1) % (fNStripZ) ;
335     Int_t col = (Int_t) TMath::Ceil((Double_t) iStrip/(fNStripZ)) -1 ;
336  
337     // Absid for 8x2-strips. Looks nice :)
338     absId = (module-1)*fNCristalsInModule +
339                   row * 2 + (col*fNCellsXInStrip + (icell - 1) / 2)*fNZ - (icell & 1 ? 1 : 0);
340  
341   }
342 */
343  
344 }
345
346 //____________________________________________________________________________
347 void AliPHOSGeoUtils::RelPosToRelId(Int_t module, Double_t x, Double_t z, Int_t * relId) const
348 {
349   //Evaluates RelId of the crystall with given coordinates
350
351   Int_t absId ;
352   RelPosToAbsId(module, x,z,absId) ;
353   AbsToRelNumbering(absId,relId) ;
354 }
355
356 //____________________________________________________________________________
357 void AliPHOSGeoUtils::RelPosInAlice(Int_t id, TVector3 & pos ) const
358 {
359   // Converts the absolute numbering into the global ALICE coordinate system
360   
361   if (!gGeoManager){
362     printf("Geo manager not initialized\n");
363     abort();
364   }
365     
366   Int_t relid[4] ;
367     
368   AbsToRelNumbering(id , relid) ;
369     
370   //construct module name
371   if(relid[1]==0){ //this is EMC
372  
373     Double_t ps[3]= {0.0,-fCryStripShift,0.}; //Position incide the crystal
374     Double_t psC[3]={0.0,0.0,0.}; //Global position
375  
376     //Shift and possibly apply misalignment corrections
377     Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/fNCellsXInStrip))*fNStripZ-
378                 (Int_t) TMath::Ceil((Double_t)relid[3]/fNCellsZInStrip) ;
379     ps[0]=((relid[2]-1)%fNCellsXInStrip-fNCellsXInStrip/2+0.5)*fCellStep ;
380     ps[2]=(-(relid[3]-1)%fNCellsZInStrip+fNCellsZInStrip/2-0.5)*fCellStep ;
381  
382     Int_t mod = relid[0] ;
383     const TGeoHMatrix * m2 = GetMatrixForStrip(mod, strip) ;
384     m2->LocalToMaster(ps,psC);
385     pos.SetXYZ(psC[0],psC[1],psC[2]) ; 
386  
387   }
388   else{
389     //first calculate position with respect to CPV plain
390     Int_t row        = relid[2] ; //offset along x axis
391     Int_t column     = relid[3] ; //offset along z axis
392     Double_t ps[3]= {0.0,fCPVBoxSizeY/2.,0.}; //Position on top of CPV
393     Double_t psC[3]={0.0,0.0,0.}; //Global position
394     pos[0] = - ( fNumberOfCPVPadsPhi/2. - row    - 0.5 ) * fPadSizePhi  ; // position of pad  with respect
395     pos[2] = - ( fNumberOfCPVPadsZ  /2. - column - 0.5 ) * fPadSizeZ  ; // of center of PHOS module
396  
397     //now apply possible shifts and rotations
398     const TGeoHMatrix *m = GetMatrixForCPV(relid[0]) ;
399     m->LocalToMaster(ps,psC);
400     pos.SetXYZ(psC[0],psC[1],-psC[2]) ; 
401   }
402
403
404 //____________________________________________________________________________
405 void AliPHOSGeoUtils::Local2Global(Int_t mod, Float_t x, Float_t z,
406                                    TVector3& globalPosition) const 
407 {
408   Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!!
409   Double_t posG[3] ;
410   const TGeoHMatrix *mPHOS = GetMatrixForModule(mod) ;
411   mPHOS->LocalToMaster(posL,posG);
412   globalPosition.SetXYZ(posG[0],posG[1],posG[2]) ;
413 }
414 //____________________________________________________________________________
415 void AliPHOSGeoUtils::Global2Local(TVector3& localPosition,
416                                    const TVector3& globalPosition,
417                                    Int_t module) const
418 {
419   // Transforms a global position to the local coordinate system
420   // of the module 
421   //Return to PHOS local system
422   Double_t posG[3]={globalPosition.X(),globalPosition.Y(),globalPosition.Z()} ;
423   Double_t posL[3]={0.,0.,0.} ;
424   const TGeoHMatrix *mPHOS = GetMatrixForModule(module) ;
425   if(mPHOS){
426     mPHOS->MasterToLocal(posG,posL);
427     localPosition.SetXYZ(posL[0],posL[1]+fCrystalShift,-posL[2]) ;  
428   }
429   else{
430     localPosition.SetXYZ(999.,999.,999.) ; //module does not exist in given configuration
431   }
432  
433 }
434 //____________________________________________________________________________
435 Bool_t AliPHOSGeoUtils::GlobalPos2RelId(TVector3 & global, Int_t * relId){
436   //Converts position in global ALICE coordinates to relId 
437   //returns false if x,z coordinates are beyond PHOS
438   //distande to PHOS surface is NOT calculated 
439   TVector3 loc ;
440   for(Int_t mod=1; mod<=fNModules; mod++){
441     Global2Local(loc,global,mod) ;
442     //If in Acceptance
443     if((TMath::Abs(loc.Z())<fXtlArrSize[2]) && (TMath::Abs(loc.X())<fXtlArrSize[0])){
444        RelPosToRelId(mod,loc.X(),loc.Z(),relId);
445        return kTRUE ;
446     }
447   }
448   return kFALSE ; 
449
450 }
451 //____________________________________________________________________________
452 Bool_t AliPHOSGeoUtils::ImpactOnEmc(const TParticle * particle,
453        Int_t & moduleNumber, Double_t & z, Double_t & x) const
454 {
455   // Tells if a particle enters PHOS and evaluates hit position
456   Double_t vtx[3]={particle->Vx(),particle->Vy(),particle->Vz()} ;
457   return ImpactOnEmc(vtx,particle->Theta(),particle->Phi(),moduleNumber,z,x);
458 }
459  
460 //____________________________________________________________________________
461 Bool_t AliPHOSGeoUtils::ImpactOnEmc(const Double_t * vtx, Double_t theta, Double_t phi, 
462                                   Int_t & moduleNumber, Double_t & z, Double_t & x) const
463 {
464   // calculates the impact coordinates on PHOS of a neutral particle  
465   // emitted in the vertex vtx[3] with direction vec(p) in the ALICE global coordinate system
466   TVector3 p(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) ;
467   return ImpactOnEmc(vtx,p,moduleNumber,z,x) ;
468
469 }
470 //____________________________________________________________________________
471 Bool_t AliPHOSGeoUtils::ImpactOnEmc(const Double_t * vtx, const TVector3 &p,
472                                   Int_t & moduleNumber, Double_t & z, Double_t & x) const
473 {
474   // calculates the impact coordinates on PHOS of a neutral particle  
475   // emitted in the vertex vtx[3] with direction theta and phi in the ALICE global coordinate system
476   TVector3 v(vtx[0],vtx[1],vtx[2]) ;
477
478   for(Int_t imod=1; imod<=fNModules ; imod++){
479     //create vector from (0,0,0) to center of crystal surface of imod module
480     Double_t tmp[3]={0.,-fCrystalShift,0.} ;
481  
482     const TGeoHMatrix *m = GetMatrixForModule(imod) ;
483     if(!m) //module does not exist in given configuration
484       continue ; 
485     Double_t posG[3]={0.,0.,0.} ;
486     m->LocalToMaster(tmp,posG);
487     TVector3 n(posG[0],posG[1],posG[2]) ; 
488     Double_t direction=n.Dot(p) ;
489     if(direction<=0.)
490       continue ; //momentum directed FROM module
491     Double_t fr = (n.Mag2()-n.Dot(v))/direction ;  
492     //Calculate direction in module plain
493     n-=v+fr*p ;
494     n*=-1. ;
495     if(TMath::Abs(TMath::Abs(n.Z())<fXtlArrSize[2]) && n.Pt()<fXtlArrSize[0]){
496       moduleNumber = imod ;
497       z=n.Z() ;
498       x=TMath::Sign(n.Pt(),n.X()) ;
499       //no need to return to local system since we calcilated distance from module center
500       //and tilts can not be significant.
501       return kTRUE ;
502     }
503   }
504   //Not in acceptance
505   x=0; z=0 ;
506   moduleNumber=0 ;
507   return kFALSE ;
508
509 }
510 //____________________________________________________________________________
511 void AliPHOSGeoUtils::GetIncidentVector(const TVector3 &vtx, Int_t module, Float_t x,Float_t z, TVector3 &vInc) const {
512   //Calculates vector pointing from vertex to current poisition in module local frame
513   //Note that PHOS local system and ALICE global have opposite z directions
514
515   Global2Local(vInc,vtx,module) ; 
516   vInc.SetXYZ(vInc.X()+x,vInc.Y(),vInc.Z()+z) ;
517 }
518 //____________________________________________________________________________
519 const TGeoHMatrix * AliPHOSGeoUtils::GetMatrixForModule(Int_t mod)const {
520   //Provides shift-rotation matrix for module mod
521
522   //If GeoManager exists, take matrixes from it
523   if(gGeoManager){
524     char path[255] ;
525     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",mod) ;
526     //    sprintf(path,"/ALIC_1/PHOS_%d",relid[0]) ;
527     if (!gGeoManager->cd(path)){
528       AliWarning(Form("Geo manager can not find path %s \n",path));
529       return 0;
530     }
531     return gGeoManager->GetCurrentMatrix();
532   }
533   if(fEMCMatrix[mod-1]){
534     return fEMCMatrix[mod-1] ;
535   }
536   else{
537     AliWarning("Can not find PHOS misalignment matrixes\n") ;
538     AliWarning("Either import TGeoManager from geometry.root or \n");
539     AliWarning("read stored matrixes from AliESD Header: \n") ;
540     AliWarning("AliPHOSGeoUtils::SetMisalMatrixes(header->GetPHOSMisalMatrix()) \n") ; 
541     return 0 ;
542   }
543   return 0 ;
544 }
545 //____________________________________________________________________________
546 const TGeoHMatrix * AliPHOSGeoUtils::GetMatrixForStrip(Int_t mod, Int_t strip)const {
547   //Provides shift-rotation matrix for strip unit of the module mod
548
549   //If GeoManager exists, take matrixes from it
550   if(gGeoManager){
551     char path[255] ;
552     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1/PSTR_%d",mod,strip) ;
553     if (!gGeoManager->cd(path)){
554       AliWarning(Form("Geo manager can not find path %s \n",path));
555       return 0 ;
556     }
557     return gGeoManager->GetCurrentMatrix();
558   } 
559   if(fStripMatrix[mod-1][strip-1]){
560     return fStripMatrix[mod-1][strip-1] ;
561   }
562   else{
563     AliWarning("Can not find PHOS misalignment matrixes\n") ;
564     AliWarning("Either import TGeoManager from geometry.root or \n");
565     AliWarning("read stored matrixes from AliESD Header: \n") ; 
566     AliWarning("AliPHOSGeoUtils::SetMisalMatrixes(header->GetPHOSMisalMatrix()) \n") ;
567     return 0 ;
568   } 
569   return 0 ;
570 }
571 //____________________________________________________________________________
572 const TGeoHMatrix * AliPHOSGeoUtils::GetMatrixForCPV(Int_t mod)const {
573   //Provides shift-rotation matrix for CPV of the module mod
574
575   //If GeoManager exists, take matrixes from it
576   if(gGeoManager){ 
577     char path[255] ;
578     //now apply possible shifts and rotations
579     sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",mod) ;
580     if (!gGeoManager->cd(path)){
581       AliWarning(Form("Geo manager can not find path %s \n",path));
582       return 0 ;
583     }
584     return gGeoManager->GetCurrentMatrix();
585   }
586   if(fCPVMatrix[mod-1]){
587     return fCPVMatrix[mod-1] ;
588   }
589   else{
590     AliWarning("Can not find PHOS misalignment matrixes\n") ;
591     AliWarning("Either import TGeoManager from geometry.root or \n");
592     AliWarning("read stored matrixes from AliESD Header: \n") ;  
593     AliWarning("AliPHOSGeoUtils::SetMisalMatrixes(header->GetPHOSMisalMatrix()) \n") ;
594     return 0 ;
595   }
596   return 0 ;
597
598 //____________________________________________________________________________
599 const TGeoHMatrix * AliPHOSGeoUtils::GetMatrixForPHOS(Int_t mod)const {
600   //Provides shift-rotation matrix for PHOS (EMC+CPV) 
601
602   //If GeoManager exists, take matrixes from it
603   if(gGeoManager){
604     char path[255] ;
605     sprintf(path,"/ALIC_1/PHOS_%d",mod) ;
606     if (!gGeoManager->cd(path)){
607       AliWarning(Form("Geo manager can not find path %s \n",path));
608       return 0 ;
609     }
610     return gGeoManager->GetCurrentMatrix();
611   }
612   if(fPHOSMatrix[mod-1]){
613     return fPHOSMatrix[mod-1] ;
614   }
615   else{
616     AliWarning("Can not find PHOS misalignment matrixes\n") ;
617     AliWarning("Either import TGeoManager from geometry.root or \n");
618     AliWarning("read stored matrixes from AliESD Header:  \n") ;   
619     AliWarning("AliPHOSGeoUtils::SetMisalMatrixes(header->GetPHOSMisalMatrix()) \n") ;
620     return 0 ;
621   }
622   return 0 ;
623 }
624 //____________________________________________________________________________
625 void AliPHOSGeoUtils::SetMisalMatrix(const TGeoHMatrix * m, Int_t mod){
626   //Fills pointers to geo matrixes
627  
628   fPHOSMatrix[mod]=m ;
629
630   //If module does not exist, make sure all its matrices are zero
631   if(m==NULL){
632     fEMCMatrix[mod]=NULL ;
633     Int_t istrip=0 ;
634     for(Int_t irow = 0; irow < fGeometryEMCA->GetNStripX(); irow ++){
635       for(Int_t icol = 0; icol < fGeometryEMCA->GetNStripZ(); icol ++){
636         fStripMatrix[mod][istrip]=NULL ;
637       }
638     } 
639     fCPVMatrix[mod]=NULL ;
640     return ;
641   }
642
643   //Calculate maxtrices for PTII
644   if(!fMisalArray)
645     fMisalArray = new TClonesArray("TGeoHMatrix",1120+10) ;
646   Int_t nr = fMisalArray->GetEntriesFast() ;
647   Double_t rotEMC[9]={1.,0.,0.,0.,0.,-1.,0.,1.,0.} ;
648   const Float_t * inthermo = fGeometryEMCA->GetInnerThermoHalfSize() ;
649   const Float_t * strip    = fGeometryEMCA->GetStripHalfSize() ;
650   const Float_t * covparams = fGeometryEMCA->GetAlCoverParams() ;
651   const Float_t * warmcov = fGeometryEMCA->GetWarmAlCoverHalfSize() ;
652   Float_t z = fGeometryCPV->GetCPVBoxSize(1) / 2. - warmcov[2] + covparams[3]-inthermo[1] ;
653   Double_t locTII[3]={0.,0.,z} ; 
654   Double_t globTII[3] ;
655
656   if (fEMCMatrix[mod] == NULL)
657     fEMCMatrix[mod] = new((*fMisalArray)[nr])TGeoHMatrix() ;
658   nr++ ;
659   fEMCMatrix[mod]->SetRotation(rotEMC) ;
660   fEMCMatrix[mod]->MultiplyLeft(fPHOSMatrix[mod]) ;
661   fPHOSMatrix[mod]->LocalToMaster(locTII,globTII) ;
662   fEMCMatrix[mod]->SetTranslation(globTII) ;
663  
664   //Now calculate ideal matrixes for strip misalignment.
665   //For the moment we can not store them in ESDHeader
666
667   Double_t loc[3]={0.,inthermo[1] - strip[1],0.} ; 
668   Double_t glob[3] ;
669
670   Int_t istrip=0 ;
671   for(Int_t irow = 0; irow < fGeometryEMCA->GetNStripX(); irow ++){
672     loc[0] = (2*irow + 1 - fGeometryEMCA->GetNStripX())* strip[0] ;
673     for(Int_t icol = 0; icol < fGeometryEMCA->GetNStripZ(); icol ++){
674       loc[2] = (2*icol + 1 - fGeometryEMCA->GetNStripZ()) * strip[2] ;
675       fEMCMatrix[mod]->LocalToMaster(loc,glob) ;
676       if (fStripMatrix[mod][istrip] == NULL)
677         fStripMatrix[mod][istrip] = new((*fMisalArray)[nr])TGeoHMatrix(*(fEMCMatrix[mod])) ; //Use same rotation as PHOS module
678       nr++ ;
679       fStripMatrix[mod][istrip]->SetTranslation(glob) ;
680       istrip++;
681     }
682   }
683  
684   //Now calculate CPV matrixes
685   const Float_t * emcParams = fGeometryEMCA->GetEMCParams() ;
686   Double_t globCPV[3] ;
687   Double_t locCPV[3]={0.,0.,- emcParams[3]} ;
688   Double_t rot[9]={1.,0.,0.,0.,0.,1.,0.,-1.,0.} ;
689
690   if (fCPVMatrix[mod] == NULL)
691     fCPVMatrix[mod] = new((*fMisalArray)[nr])TGeoHMatrix() ;
692   nr++ ;
693   fCPVMatrix[mod]->SetRotation(rot) ;
694   fCPVMatrix[mod]->MultiplyLeft(fPHOSMatrix[mod]) ;
695   fCPVMatrix[mod]->ReflectY(kFALSE) ;
696   fPHOSMatrix[mod]->LocalToMaster(locCPV,globCPV) ;
697   fCPVMatrix[mod]->SetTranslation(globCPV) ;
698
699 }
700