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