set minimum pT for dicing (M. Verweij)
[u/mrichter/AliRoot.git] / JETAN / AliJetDummyGeo.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 //
17 // Temporarily added to define part of the EMCal geometry
18 // necessary for the jet finder
19 //
20 // M. Estienne
21 // Magali.Estienne@cern.ch
22 //
23
24 #include <Riostream.h>
25 #include <assert.h>
26 #include <TList.h>
27
28 // --- Root header files ---
29 #include <TVector3.h>
30 #include <TGeoManager.h>
31 #include <TGeoMatrix.h>
32 #include <TGeoNode.h>
33
34 // --- AliRoot header files ---
35 #include "AliLog.h"
36 #include "AliJetDummyGeo.h"
37 #include "AliJetDummyShishKebabTrd1Module.h"
38
39 using std::cout;
40 using std::endl;
41 ClassImp(AliJetDummyGeo)
42
43 AliJetDummyGeo::AliJetDummyGeo():
44   fArm1EtaMin(-0.7), 
45   fArm1EtaMax(+0.7), 
46   fArm1PhiMin(80.0), 
47   fArm1PhiMax(200.0),
48   fNumberOfSuperModules(12), 
49   fSteelFrontThick(0.0), 
50   fLateralSteelStrip(0.01),
51   fIPDistance(428.0), 
52   fPhiGapForSM(2.), 
53   fNPhi(12),       
54   fNZ(24),      
55   fPhiModuleSize(12.26 - fPhiGapForSM / Float_t(fNPhi)), 
56   fEtaModuleSize(fPhiModuleSize),
57   fNPHIdiv(2), 
58   fNETAdiv(2), 
59   fPhiTileSize(fPhiModuleSize/Double_t(fNPHIdiv) - fLateralSteelStrip),
60   fEtaTileSize(fEtaModuleSize/Double_t(fNETAdiv) - fLateralSteelStrip),
61   fNECLayers(77), 
62   fECScintThick(0.16), 
63   fECPbRadThickness(0.16), 
64   fSampling(12.327),
65   fTrd1Angle(1.5), 
66   fNCellsInModule(fNPHIdiv*fNETAdiv), 
67   fNCellsInSupMod(fNCellsInModule*fNPhi*fNZ), 
68   fNCells(fNCellsInSupMod*fNumberOfSuperModules-fNCellsInSupMod), 
69   fLongModuleSize(fNECLayers*(fECScintThick + fECPbRadThickness)), 
70   f2Trd1Dx2(fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.)), 
71   fShellThickness(TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2)+fSteelFrontThick),
72   fZLength(2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax)),
73   fEtaMaxOfTRD1(0.), 
74   fPhiBoundariesOfSM(0),
75   fPhiCentersOfSM(0),
76   fCentersOfCellsEtaDir(0),
77   fCentersOfCellsXDir(0),
78   fCentersOfCellsPhiDir(0),
79   fEtaCentersOfCells(0),
80   fPhiCentersOfCells(0),
81   fShishKebabTrd1Modules(0),
82   fDebug(0)
83 {
84   // Constructor
85
86   // Local coordinates
87   fParSM[0] = GetShellThickness()/2.;        
88   fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
89   fParSM[2] = 350./2.;
90
91   fEnvelop[0] = fIPDistance; // mother volume inner radius
92   fEnvelop[1] = fIPDistance + fShellThickness; // mother volume outer r.
93   fEnvelop[2] = 1.00001*fZLength; // add some padding for mother volume. 
94
95   // SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006 
96   fPhiBoundariesOfSM.Set(fNumberOfSuperModules);
97   fPhiCentersOfSM.Set(fNumberOfSuperModules/2);
98   fPhiBoundariesOfSM[0] = TMath::PiOver2() - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules)
99   fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance);
100   fPhiCentersOfSM[0]    = TMath::PiOver2();
101   for(Int_t i=1; i<=4; i++) { // from 2th ro 9th
102     fPhiBoundariesOfSM[2*i]   = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i;
103     fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i;
104     fPhiCentersOfSM[i]        = fPhiCentersOfSM[0]    + 20.*TMath::DegToRad()*i;
105   }
106   fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad();
107   fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance);
108   fPhiCentersOfSM[5]     = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.; 
109
110   Int_t nphism  = GetNumberOfSuperModules()/2;
111   Double_t dphi = (GetArm1PhiMax() - GetArm1PhiMin())/nphism;
112   Double_t rpos = (GetEnvelop(0) + GetEnvelop(1))/2.;
113   Double_t phi, phiRad, xpos, ypos, zpos;
114   for(Int_t i=0; i<nphism; i++){
115     phi    = GetArm1PhiMin() + dphi*(2*i+1)/2.; // phi= 90, 110, 130, 150, 170, 190
116     phiRad = phi*TMath::Pi()/180.;
117     xpos = rpos * TMath::Cos(phiRad);
118     ypos = rpos * TMath::Sin(phiRad);
119     zpos = fParSM[2];
120     if(i==5) {
121       xpos += (fParSM[1]/2. * TMath::Sin(phiRad)); 
122       ypos -= (fParSM[1]/2. * TMath::Cos(phiRad));
123     }
124     // pozitive z
125     Int_t ind = 2*i;
126     TGeoRotation *geoRot0 = new TGeoRotation("geoRot0", 90.0, phi, 90.0, 90.0+phi, 0.0, 0.0);
127     fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind),
128                                           xpos,ypos, zpos, geoRot0);
129     // negaive z
130     ind++;
131     Double_t phiy = 90. + phi + 180.;
132     if(phiy>=360.) phiy -= 360.;
133     TGeoRotation *geoRot1 = new TGeoRotation("geoRot1", 90.0, phi, 90.0, phiy, 180.0, 0.0);
134     fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind),
135                                           xpos,ypos,-zpos, geoRot1);
136   } // for
137
138   CreateListOfTrd1Modules();
139
140   if(fDebug > 0){
141     for(Int_t i=0; i<6; i++){
142       cout << "fMatrixOfSM[" << i << "]: " << fMatrixOfSM[i] << endl;
143     }
144     cout << "fArm1EtaMin: " << fArm1EtaMin << endl;
145     cout << "fArm1EtaMax: " << fArm1EtaMax << endl;
146     cout << "fArm1PhiMin: " << fArm1PhiMin << endl;
147     cout << "fArm1PhiMax: " << fArm1PhiMax << endl;
148     cout << "fNumberOfSuperModules: " << fNumberOfSuperModules << endl;
149     cout << "fSteelFrontThick: " << fSteelFrontThick << endl;
150     cout << "fIPDistance: " << fIPDistance  << endl;
151     cout << "fZLength: " << fZLength << endl;
152     cout << "fPhiGapForSM: " << fPhiGapForSM << endl;
153     cout << "fNPhi: " << fNPhi << endl;
154     cout << "fNZ: " << fNZ << endl;
155     cout << "fPhiModuleSize: " << fPhiModuleSize << endl;
156     cout << "fEtaModuleSize: " << fEtaModuleSize << endl;
157     cout << "fNPHIdiv: " << fNPHIdiv << endl;
158     cout << "fNETAdiv: " << fNETAdiv << endl;
159     cout << "fNECLayers: " << fNECLayers << endl;
160     cout << "fECScintThick: " <<  fECScintThick << endl;
161     cout << "fECPbRadThickness: " << fECPbRadThickness << endl;
162     cout << "fSampling: " << fSampling << endl;
163     cout << "fTrd1Angle: " << fTrd1Angle << endl;
164     cout << "fNCellsInModule: " << fNCellsInModule << endl;
165     cout << "fNCellsInSupMod: " << fNCellsInSupMod << endl;
166     cout << "fNCells: " << fNCells << endl;
167     cout << "fLongModuleSize: " <<  fLongModuleSize << endl;
168     cout << "f2Trd1Dx2: " << f2Trd1Dx2 << endl;
169     cout << "fShellThickness: " << fShellThickness << endl;
170     cout << "fEtaMaxOfTRD1: " << fEtaMaxOfTRD1 << endl;
171   }
172 }
173
174 AliJetDummyGeo::AliJetDummyGeo(const AliJetDummyGeo& geom):
175   TObject(geom),  
176   fArm1EtaMin(geom.fArm1EtaMin), 
177   fArm1EtaMax(geom.fArm1EtaMax), 
178   fArm1PhiMin(geom.fArm1PhiMin), 
179   fArm1PhiMax(geom.fArm1PhiMax),
180   fNumberOfSuperModules(geom.fNumberOfSuperModules), 
181   fSteelFrontThick(geom.fSteelFrontThick), 
182   fLateralSteelStrip(geom.fLateralSteelStrip),
183   fIPDistance(geom.fIPDistance), 
184   fPhiGapForSM(geom.fPhiGapForSM), 
185   fNPhi(geom.fNPhi),       
186   fNZ(geom.fNZ),      
187   fPhiModuleSize(geom.fPhiModuleSize), 
188   fEtaModuleSize(geom.fEtaModuleSize),
189   fNPHIdiv(geom.fNPHIdiv), 
190   fNETAdiv(geom.fNETAdiv), 
191   fPhiTileSize(geom.fPhiTileSize),
192   fEtaTileSize(geom.fEtaTileSize),
193   fNECLayers(geom.fNECLayers), 
194   fECScintThick(geom.fECScintThick), 
195   fECPbRadThickness(geom.fECPbRadThickness), 
196   fSampling(geom.fSampling),
197   fTrd1Angle(geom.fTrd1Angle), 
198   fNCellsInModule(geom.fNCellsInModule), 
199   fNCellsInSupMod(geom.fNCellsInSupMod), 
200   fNCells(geom.fNCells),
201   fLongModuleSize(geom.fLongModuleSize),
202   f2Trd1Dx2(geom.f2Trd1Dx2),
203   fShellThickness(geom.fShellThickness),
204   fZLength(geom.fZLength),
205   fEtaMaxOfTRD1(geom.fEtaMaxOfTRD1),
206   fPhiBoundariesOfSM(geom.fPhiBoundariesOfSM),
207   fPhiCentersOfSM(geom.fPhiCentersOfSM),
208   fCentersOfCellsEtaDir(geom.fCentersOfCellsEtaDir),
209   fCentersOfCellsXDir(geom.fCentersOfCellsXDir),
210   fCentersOfCellsPhiDir(geom.fCentersOfCellsPhiDir),
211   fEtaCentersOfCells(geom.fEtaCentersOfCells),
212   fPhiCentersOfCells(geom.fPhiCentersOfCells),
213   fShishKebabTrd1Modules(geom.fShishKebabTrd1Modules),
214   fDebug(geom.fDebug)
215 {
216   // Constructor
217   // Local coordinates
218   fParSM[0] = GetShellThickness()/2.;        
219   fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
220   fParSM[2] = 350./2.;
221
222   fEnvelop[0] = fIPDistance; // mother volume inner radius
223   fEnvelop[1] = fIPDistance + fShellThickness; // mother volume outer r.
224   fEnvelop[2] = 1.00001*fZLength; // add some padding for mother volume. 
225
226   // SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006 
227   fPhiBoundariesOfSM.Set(fNumberOfSuperModules);
228   fPhiCentersOfSM.Set(fNumberOfSuperModules/2);
229   fPhiBoundariesOfSM[0] = TMath::PiOver2() - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules)
230   fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance);
231   fPhiCentersOfSM[0]    = TMath::PiOver2();
232   for(Int_t i=1; i<=4; i++) { // from 2th ro 9th
233     fPhiBoundariesOfSM[2*i]   = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i;
234     fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i;
235     fPhiCentersOfSM[i]        = fPhiCentersOfSM[0]    + 20.*TMath::DegToRad()*i;
236   }
237   fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad();
238   fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance);
239   fPhiCentersOfSM[5]     = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.; 
240
241   Int_t nphism  = GetNumberOfSuperModules()/2;
242   Double_t dphi = (GetArm1PhiMax() - GetArm1PhiMin())/nphism;
243   Double_t rpos = (GetEnvelop(0) + GetEnvelop(1))/2.;
244   Double_t phi, phiRad, xpos, ypos, zpos;
245   for(Int_t i=0; i<nphism; i++){
246     phi    = GetArm1PhiMin() + dphi*(2*i+1)/2.; // phi= 90, 110, 130, 150, 170, 190
247     phiRad = phi*TMath::Pi()/180.;
248     xpos = rpos * TMath::Cos(phiRad);
249     ypos = rpos * TMath::Sin(phiRad);
250     zpos = fParSM[2];
251     if(i==5) {
252       xpos += (fParSM[1]/2. * TMath::Sin(phiRad)); 
253       ypos -= (fParSM[1]/2. * TMath::Cos(phiRad));
254     }
255     // pozitive z
256     Int_t ind = 2*i;
257     TGeoRotation *geoRot0 = new TGeoRotation("geoRot0", 90.0, phi, 90.0, 90.0+phi, 0.0, 0.0);
258     fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind),
259                                           xpos,ypos, zpos, geoRot0);
260     // negaive z
261     ind++;
262     Double_t phiy = 90. + phi + 180.;
263     if(phiy>=360.) phiy -= 360.;
264     TGeoRotation *geoRot1 = new TGeoRotation("geoRot1", 90.0, phi, 90.0, phiy, 180.0, 0.0);
265     fMatrixOfSM[ind] = new TGeoCombiTrans(Form("EmcalSM%2.2i",ind),
266                                           xpos,ypos,-zpos, geoRot1);
267
268     delete geoRot0;
269     delete geoRot1;
270
271   } // for
272
273   CreateListOfTrd1Modules();
274
275 }
276
277 //------------------------------------------------------------------------------------
278 AliJetDummyGeo::~AliJetDummyGeo()
279 {
280   // Destructor
281   // delete [] fMatrixOfSM;
282 }
283
284 //------------------------------------------------------------------------------------
285 void AliJetDummyGeo::EtaPhiFromIndex(Int_t absId, Float_t& eta, Float_t& phi) const
286 {
287   // Nov 16, 2006- float to double
288   // version for TRD1 only
289   static TVector3 vglob;
290   GetGlobal(absId, vglob);
291   eta = vglob.Eta();
292   phi = vglob.Phi();
293
294 }
295
296 //------------------------------------------------------------------------------------
297 void AliJetDummyGeo::GetGlobal(const Double_t *loc, Double_t *glob, Int_t ind) const
298 {
299   // Figure out the global numbering of a given supermodule from the
300   // local numbering
301   // Alice numbering - Jun 03,2006
302   //  if(fMatrixOfSM[0] == 0) GetTransformationForSM();
303
304   if(ind>=0 && ind < GetNumberOfSuperModules()) {
305     fMatrixOfSM[ind]->LocalToMaster(loc, glob);
306   }
307 }
308
309 //------------------------------------------------------------------------------------
310 void AliJetDummyGeo::GetGlobal(Int_t absId , Double_t glob[3]) const
311
312   // Alice numbering scheme - Jun 03, 2006
313   static Int_t nSupMod, nModule, nIphi, nIeta;
314   static Double_t loc[3];
315
316   glob[0]=glob[1]=glob[2]=0.0; // bad case
317   if(RelPosCellInSModule(absId, loc)) {
318     GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
319     fMatrixOfSM[nSupMod]->LocalToMaster(loc, glob);
320   }
321 }
322
323 //------------------------------------------------------------------------------------
324 void AliJetDummyGeo::GetGlobal(Int_t absId , TVector3 &vglob) const
325
326   // Alice numbering scheme - Jun 03, 2006
327   static Double_t glob[3];
328
329   GetGlobal(absId, glob);
330   vglob.SetXYZ(glob[0], glob[1], glob[2]);
331
332 }
333
334 //------------------------------------------------------------------------------------
335 Bool_t AliJetDummyGeo::RelPosCellInSModule(Int_t absId, Double_t loc[3]) const
336 {
337   // Alice numbering scheme - Jun 03, 2006
338   loc[0] = loc[1] = loc[2]=0.0;
339   if(RelPosCellInSModule(absId, loc[0],loc[1],loc[2])) {
340     return kTRUE;
341   }
342   return kFALSE;
343 }
344
345 //------------------------------------------------------------------------------------
346 Bool_t AliJetDummyGeo::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const
347 {
348   // Look to see what the relative position inside a given cell is
349   // for a recpoint.
350   // Alice numbering scheme - Jun 08, 2006
351   // In:
352   // absId   - cell is as in Geant,     0<= absId   < fNCells;
353   // OUT:
354   // xr,yr,zr - x,y,z coordinates of cell with absId inside SM 
355
356   // Shift index taking into account the difference between standard SM 
357   // and SM of half size in phi direction
358   const Int_t phiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
359   static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
360   if(!CheckAbsCellId(absId)) return kFALSE;
361
362   GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
363   GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta); 
364  
365   xr = fCentersOfCellsXDir.At(ieta);
366   zr = fCentersOfCellsEtaDir.At(ieta);
367
368   if(nSupMod<10) {
369     yr = fCentersOfCellsPhiDir.At(iphi);
370   } else {
371     yr = fCentersOfCellsPhiDir.At(iphi + phiIndexShift);
372   }
373
374   return kTRUE;
375 }
376
377 //------------------------------------------------------------------------------------
378 Bool_t  AliJetDummyGeo::CheckAbsCellId(Int_t absId) const
379
380   // May 31, 2006; only trd1 now
381   if(absId<0 || absId >= fNCells) return kFALSE;
382   else                            return kTRUE;
383 }
384
385 //------------------------------------------------------------------------------------
386 Bool_t AliJetDummyGeo::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nModule,Int_t &nIphi,Int_t &nIeta) const
387
388   // 21-sep-04; 19-oct-05;
389   // May 31, 2006; ALICE numbering scheme:
390   // 
391   // In:
392   // absId   - cell is as in Geant,     0<= absId   < fNCells;
393   // Out:
394   // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
395   // nModule  - module number in SM,    0<= nModule  < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
396   // nIphi   - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv; 
397   // nIeta   - cell number in eta driection inside module; 0<= nIeta < fNETAdiv; 
398   // 
399   static Int_t tmp=0, sm10=0;
400   if(!CheckAbsCellId(absId)) return kFALSE;
401
402   sm10 = fNCellsInSupMod*10;
403   if(absId >= sm10) { // 110 degree case; last two supermodules  
404     nSupMod = (absId-sm10) / (fNCellsInSupMod/2) + 10;
405     tmp     = (absId-sm10) % (fNCellsInSupMod/2);
406   } else {
407     nSupMod = absId / fNCellsInSupMod;
408     tmp     = absId % fNCellsInSupMod;
409   }
410
411   nModule  = tmp / fNCellsInModule;
412   tmp     = tmp % fNCellsInModule;
413   nIphi   = tmp / fNPHIdiv;
414   nIeta   = tmp % fNPHIdiv;
415
416   return kTRUE;
417 }
418
419 //------------------------------------------------------------------------------------
420 void AliJetDummyGeo::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta, Int_t &iphi, Int_t &ieta) const
421
422   // 
423   // Added nSupMod; Nov 25, 05
424   // Alice numbering scheme  - Jun 01,2006 
425   // IN:
426   // nSupMod - super module(SM) number, 0<= nSupMod < fNumberOfSuperModules;
427   // nModule  - module number in SM,     0<= nModule  < fNCellsInSupMod/fNCellsInSupMod or(/2) for tow last SM (10th and 11th);
428   // nIphi   - cell number in phi driection inside module; 0<= nIphi < fNPHIdiv; 
429   // nIeta   - cell number in eta driection inside module; 0<= nIeta < fNETAdiv; 
430   // 
431  // OUT:
432   // ieta, iphi - indexes of cell(tower) in two dimensional grid of SM
433   // ieta - have to change from 0 to (fNZ*fNETAdiv-1)
434   // iphi - have to change from 0 to (fNPhi*fNPHIdiv-1 or fNPhi*fNPHIdiv/2-1)
435   //
436   static Int_t iphim, ietam;
437
438   GetModulePhiEtaIndexInSModule(nSupMod,nModule, iphim, ietam); 
439   //  ieta  = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM) 
440   ieta  = ietam*fNETAdiv + (fNETAdiv - 1 - nIeta); // x(module) = -z(SM) 
441   iphi  = iphim*fNPHIdiv + nIphi;     // y(module) =  y(SM) 
442
443 }
444
445
446 //------------------------------------------------------------------------------------
447 void AliJetDummyGeo::GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule,  Int_t &iphim, Int_t &ietam) const
448
449   // added nSupMod; - 19-oct-05 !
450   // Alice numbering scheme        - Jun 01,2006 
451   // ietam, iphi - indexes of module in two dimensional grid of SM
452   // ietam - have to change from 0 to fNZ-1
453   // iphim - have to change from 0 to nphi-1 (fNPhi-1 or fNPhi/2-1)
454   static Int_t nphi;
455
456   if(nSupMod>=10) nphi = fNPhi/2;
457   else            nphi = fNPhi;
458
459   ietam = nModule/nphi;
460   iphim = nModule%nphi;
461 }
462
463 //------------------------------------------------------------------------------------
464 Bool_t AliJetDummyGeo::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi, Int_t &absId) const
465 {
466   // Nov 17,2006
467   // stay here - phi problem as usual 
468   static Int_t nSupMod, i, ieta, iphi, etaShift, nphi;
469   static Double_t absEta=0.0, d=0.0, dmin=0.0, phiLoc;
470   absId = nSupMod = - 1;
471   if(SuperModuleNumberFromEtaPhi(eta, phi, nSupMod)) {
472     // phi index first
473     phi    = TVector2::Phi_0_2pi(phi);
474     phiLoc = phi - fPhiCentersOfSM[nSupMod/2];
475     nphi   = fPhiCentersOfCells.GetSize();
476     if(nSupMod>=10) {
477       phiLoc = phi - 190.*TMath::DegToRad();
478       nphi /= 2;
479     }
480
481     dmin   = TMath::Abs(fPhiCentersOfCells[0]-phiLoc);
482     iphi   = 0;
483     for(i=1; i<nphi; i++) {
484       d = TMath::Abs(fPhiCentersOfCells[i] - phiLoc);
485       if(d < dmin) {
486         dmin = d;
487         iphi = i;
488       }
489     }
490     // odd SM are turned with respect of even SM - reverse indexes
491
492     // eta index
493     absEta   = TMath::Abs(eta);
494     etaShift = iphi*fCentersOfCellsEtaDir.GetSize();
495     dmin     = TMath::Abs(fEtaCentersOfCells[etaShift]-absEta);
496     ieta     = 0;
497     for(i=1; i<fCentersOfCellsEtaDir.GetSize(); i++) {
498       d = TMath::Abs(fEtaCentersOfCells[i+etaShift] - absEta);
499       if(d < dmin) {
500         dmin = d;
501         ieta = i;
502       }
503     }
504
505     if(eta<0) iphi = (nphi-1) - iphi;
506     absId = GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);
507
508     return kTRUE;
509   }
510   return kFALSE;
511 }
512
513 //------------------------------------------------------------------------------------
514 Bool_t AliJetDummyGeo::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const
515
516   // Return false if phi belongs a phi cracks between SM
517  
518   static Int_t i;
519
520   if(TMath::Abs(eta) > fEtaMaxOfTRD1) return kFALSE;
521
522   phi = TVector2::Phi_0_2pi(phi); // move phi to (0,2pi) boundaries
523   for(i=0; i<6; i++) {
524     if(phi>=fPhiBoundariesOfSM[2*i] && phi<=fPhiBoundariesOfSM[2*i+1]) {
525       nSupMod = 2*i;
526       if(eta < 0.0) nSupMod++;
527       return kTRUE;
528     }
529   }
530   return kFALSE;
531 }
532
533 //------------------------------------------------------------------------------------
534 Int_t  AliJetDummyGeo::GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const
535 {
536   // Transition from super module number(nSupMod) and cell indexes (ieta,iphi) to absId
537   static Int_t ietam, iphim, nModule;
538   static Int_t nIeta, nIphi; // cell indexes in module
539
540   GetModuleIndexesFromCellIndexesInSModule(nSupMod, iphi, ieta, ietam, iphim, nModule);
541
542   nIeta = ieta%fNETAdiv;
543   nIeta = fNETAdiv - 1 - nIeta;
544   nIphi = iphi%fNPHIdiv;
545
546   return GetAbsCellId(nSupMod, nModule, nIphi, nIeta);
547 }
548
549 //------------------------------------------------------------------------------------
550 void  AliJetDummyGeo::GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta, 
551                         Int_t &iphim, Int_t &ietam, Int_t &nModule) const
552 {
553   // Transition from cell indexes (ieta,iphi) to module indexes (ietam,iphim, nModule)
554   static Int_t nphi;
555   nphi  = GetNumberOfModuleInPhiDirection(nSupMod);  
556
557   ietam  = ieta/fNETAdiv;
558   iphim  = iphi/fNPHIdiv;
559   nModule = ietam * nphi + iphim; 
560 }
561
562 //------------------------------------------------------------------------------------
563 Int_t AliJetDummyGeo::GetAbsCellId(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta) const
564
565   // 27-aug-04; 
566   // corr. 21-sep-04; 
567   //       13-oct-05; 110 degree case
568   // May 31, 2006; ALICE numbering scheme:
569   // 0 <= nSupMod < fNumberOfSuperModules
570   // 0 <= nModule  < fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1)
571   // 0 <= nIphi   < fNPHIdiv
572   // 0 <= nIeta   < fNETAdiv
573   // 0 <= absid   < fNCells
574   static Int_t id=0; // have to change from 0 to fNCells-1
575   if(nSupMod >= 10) { // 110 degree case; last two supermodules
576     id  = fNCellsInSupMod*10 + (fNCellsInSupMod/2)*(nSupMod-10);
577   } else {
578     id  = fNCellsInSupMod*nSupMod;
579   }
580   id += fNCellsInModule *nModule;
581   id += fNPHIdiv *nIphi;
582   id += nIeta;
583   if(id<0 || id >= fNCells) {
584     id = -TMath::Abs(id); // if negative something wrong
585   }
586   return id;
587 }
588
589 //------------------------------------------------------------------------------------
590 Bool_t AliJetDummyGeo::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const
591 {
592   // 0<= nPhiSec <=4; phi in rad
593   // 0;  gap boundaries between  0th&2th  | 1th&3th SM
594   // 1;  gap boundaries between  2th&4th  | 3th&5th SM
595   // 2;  gap boundaries between  4th&6th  | 5th&7th SM
596   // 3;  gap boundaries between  6th&8th  | 7th&9th SM
597   // 4;  gap boundaries between  8th&10th | 9th&11th SM
598   if(nPhiSec<0 || nPhiSec >4) return kFALSE; 
599   phiMin = fPhiBoundariesOfSM[2*nPhiSec+1];
600   phiMax = fPhiBoundariesOfSM[2*nPhiSec+2];
601   return kTRUE; 
602 }
603
604 //------------------------------------------------------------------------------------
605 void AliJetDummyGeo::CreateListOfTrd1Modules()
606 {
607   // Generate the list of Trd1 modules
608   // which will make up the EMCAL
609   // geometry
610     printf("CreateListOfTrd1Modules() \n");
611     
612   AliJetDummyShishKebabTrd1Module *mod=0, *mTmp=0; // current module
613   if(fShishKebabTrd1Modules == 0) {
614     fShishKebabTrd1Modules = new TList;
615     fShishKebabTrd1Modules->SetName("ListOfTRD1");
616     for(Int_t iz=0; iz< GetNZ(); iz++) { 
617       if(iz==0) { 
618         mod  = new AliJetDummyShishKebabTrd1Module(TMath::Pi()/2.,this);
619       } else {
620         mTmp  = new AliJetDummyShishKebabTrd1Module(*mod);
621         mod   = mTmp;
622       }
623       fShishKebabTrd1Modules->Add(mod);
624     }
625   } else {
626     printf(" Already exits : ");
627   }
628   mod = (AliJetDummyShishKebabTrd1Module*)fShishKebabTrd1Modules->At(fShishKebabTrd1Modules->GetSize()-1);
629   fEtaMaxOfTRD1 = mod->GetMaxEtaOfModule(0);
630
631   //  printf(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n", 
632   //              fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1);
633
634   // Feb 20,2006;
635   // Jun 01, 2006 - ALICE numbering scheme
636   // define grid for cells in eta(z) and x directions in local coordinates system of SM
637   // Works just for 2x2 case only -- ?? start here
638   // 
639   //
640   // Define grid for cells in phi(y) direction in local coordinates system of SM
641   // as for 2X2 as for 3X3 - Nov 8,2006
642   // 
643
644   //  printf(" Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize());
645
646   Int_t ind=0; // this is phi index
647   Int_t ieta=0, nModule=0, iphiTemp;
648   Double_t xr, zr, theta, phi, eta, r, x,y;
649   xr = 0.;
650   zr = 0.;
651   
652   TVector3 vglob;
653   Double_t ytCenterModule=0.0, ytCenterCell=0.0;
654
655   fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
656   fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
657
658   Double_t R0 = GetIPDistance() + GetLongModuleSize()/2.;
659   for(Int_t it=0; it<fNPhi; it++) { // cycle on modules
660     ytCenterModule = -fParSM[1] + fPhiModuleSize*(2*it+1)/2;  // center of module
661     for(Int_t ic=0; ic<fNPHIdiv; ic++) { // cycle on cells in module
662       if(fNPHIdiv==2) {
663         ytCenterCell = ytCenterModule + fPhiTileSize *(2*ic-1)/2.;
664       } else if(fNPHIdiv==3){
665         ytCenterCell = ytCenterModule + fPhiTileSize *(ic-1);
666       } else if(fNPHIdiv==1){
667         ytCenterCell = ytCenterModule;
668       }
669       fCentersOfCellsPhiDir.AddAt(ytCenterCell,ind);
670       // Define grid on phi direction
671       // Grid is not the same for different eta bin;
672       // Effect is small but is still here
673       phi = TMath::ATan2(ytCenterCell, R0);
674       fPhiCentersOfCells.AddAt(phi, ind);
675
676       //      printf(" ind %2.2i : y %8.3f ", ind, fCentersOfCellsPhiDir.At(ind)); 
677       ind++;
678     }
679   }
680
681   fCentersOfCellsEtaDir.Set(fNZ *fNETAdiv);
682   fCentersOfCellsXDir.Set(fNZ *fNETAdiv);
683   fEtaCentersOfCells.Set(fNZ *fNETAdiv * fNPhi*fNPHIdiv);
684   if(fDebug>1) AliInfo(Form(" Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize()));
685   for(Int_t it=0; it<fNZ; it++) {
686     AliJetDummyShishKebabTrd1Module *trd1 = GetShishKebabModule(it);
687     nModule = fNPhi*it;
688     for(Int_t ic=0; ic<fNETAdiv; ic++) {
689       if(fNPHIdiv==2) {
690         trd1->GetCenterOfCellInLocalCoordinateofSM(ic, xr, zr);      // case of 2X2
691         GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta); 
692       } if(fNPHIdiv==3) {
693         trd1->GetCenterOfCellInLocalCoordinateofSM3X3(ic, xr, zr);  // case of 3X3
694         GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta); 
695       } if(fNPHIdiv==1) {
696         trd1->GetCenterOfCellInLocalCoordinateofSM1X1(xr, zr);      // case of 1X1
697         GetCellPhiEtaIndexInSModule(0, nModule, 0, ic, iphiTemp, ieta); 
698       }
699       fCentersOfCellsXDir.AddAt(Float_t(xr) - fParSM[0],ieta);
700       fCentersOfCellsEtaDir.AddAt(Float_t(zr) - fParSM[2],ieta);
701       // Define grid on eta direction for each bin in phi
702       for(Int_t iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
703         x = xr + trd1->GetRadius();
704         y = fCentersOfCellsPhiDir[iphi];
705         r = TMath::Sqrt(x*x + y*y + zr*zr);
706         theta = TMath::ACos(zr/r);
707         eta   = AliJetDummyShishKebabTrd1Module::ThetaToEta(theta);
708         //        ind   = ieta*fCentersOfCellsPhiDir.GetSize() + iphi;
709         ind   = iphi*fCentersOfCellsEtaDir.GetSize() + ieta;
710         fEtaCentersOfCells.AddAt(eta, ind);
711       }
712       //printf(" ieta %i : xr + trd1->GetRadius() %f : zr %f : eta %f \n", ieta, xr + trd1->GetRadius(), zr, eta);
713     }
714   }
715
716   if(fDebug>1){
717     for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
718       AliInfo(Form(" ind %2.2i : z %8.3f : x %8.3f", i+1, 
719                    fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i)));
720     }
721   }
722 }
723
724 //------------------------------------------------------------------------------------
725 AliJetDummyShishKebabTrd1Module* AliJetDummyGeo::GetShishKebabModule(Int_t neta)
726 {
727   // This method was too long to be included in the header file - the
728   // rule checker complained about it's length, so we move it here. It
729   // returns the shishkebabmodule at a given eta index point.
730
731   static AliJetDummyShishKebabTrd1Module* trd1=0;
732   if(fShishKebabTrd1Modules && neta>=0 && neta<fShishKebabTrd1Modules->GetSize()) {
733     trd1 = (AliJetDummyShishKebabTrd1Module*)fShishKebabTrd1Modules->At(neta);
734   } else trd1 = 0;
735   return trd1;
736 }
737
738 //------------------------------------------------------------------------------------
739 void  AliJetDummyGeo::GetTransformationForSM()
740 {
741   // Uses the geometry manager to load the transformation matrix
742   // for the supermodules
743   // Unused after 19 Jan, 2007 - keep for compatibility; 
744
745   return;
746   static Bool_t transInit=kFALSE;
747   if(transInit) return;
748
749   Int_t i=0;
750   if(gGeoManager == 0) {
751     Info("CreateTransformationForSM() "," Load geometry : TGeoManager::Import()");
752     assert(0);
753   }
754
755   TGeoNode *tn = gGeoManager->GetTopNode();
756   TGeoNode *node=0, *xen1 = 0;
757   for(i=0; i<tn->GetNdaughters(); i++) {
758     node = tn->GetDaughter(i);
759     TString ns(node->GetName());
760     if(ns.Contains(GetNameOfEMCALEnvelope())) {
761       xen1 = node;
762       break;
763     }
764   }
765
766   if(!xen1) {
767     Info("CreateTransformationForSM() "," geometry has not EMCAL envelope with name %s", 
768     GetNameOfEMCALEnvelope());
769     assert(0);
770   }
771   AliInfo(Form(" i %i : EMCAL Envelope is %s : #SM %i \n", i, xen1->GetName(), xen1->GetNdaughters()));
772   for(i=0; i<xen1->GetNdaughters(); i++) {
773     TGeoNodeMatrix *sm = (TGeoNodeMatrix*)xen1->GetDaughter(i);
774     fMatrixOfSM[i] = sm->GetMatrix();
775   }
776   AliInfo(Form("transInit %d: ", transInit));
777   transInit = kTRUE;
778 }
779