Avoid message E-TClonesArray::At: during digitization due to try to access non existi...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALGeometry.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$*/
17
18 //_________________________________________________________________________
19 // Geometry class  for EMCAL : singleton  
20 // EMCAL consists of layers of scintillator and lead
21 // with scintillator fiber arranged as "shish-kebab" skewers 
22 // Places the the Barrel Geometry of The EMCAL at Midrapidity
23 // between 80 and 180(or 190) degrees of Phi and
24 // -0.7 to 0.7 in eta 
25 //
26 //     EMCAL geometry tree:
27 //     EMCAL -> superModule -> module -> tower(cell)
28 //     Indexes
29 //     absId -> nSupMod     -> nModule -> (nIphi,nIeta)
30 //
31 //   Name choices: 
32 //   EMCAL_PDC06 (geometry used for PDC06 simulations, kept for backward compatibility)
33 //      = equivalent to SHISH_77_TRD1_2X2_FINAL_110DEG in old notation
34 //   EMCAL_COMPLETE (geometry for expected complete detector)
35 //      = equivalent to SHISH_77_TRD1_2X2_FINAL_110DEG scTh=0.176 pbTh=0.144
36 //          in old notation
37 //   EMCAL_FIRSTYEARV1 - geometry for December 2009 to December 2010 run period; 
38 //                fixed bug for positions of modules inside SM
39 //                (first module has tilt 0.75 degree);
40 //                the sizes updated with last information from production
41 //                drawing (end of October 2010). 
42 //      
43 //   EMCAL_COMPLETEV1: Same fixes as FIRSTYEAR and 10 SM instead of 10+2 half SM
44 //
45 //   EMCAL_WSUC (Wayne State test stand)
46 //      = no definite equivalent in old notation, was only used by
47 //          Aleksei, but kept for testing purposes
48 //
49 //   etc.
50 //
51 //
52 //
53 //*-- Author: Sahal Yacoob (LBL / UCT)
54 //     and  : Yves Schutz (SUBATECH)
55 //     and  : Jennifer Klay (LBL)
56 //     and  : Alexei Pavlinov (WSU) 
57 //
58
59 //--- Root header files ---
60 #include <TVector2.h>
61 #include <TVector3.h>
62 //-- ALICE Headers.
63 #include "AliLog.h"
64
65 // // --- EMCAL headers
66 #include "AliEMCALGeometry.h"
67 #include "AliEMCALShishKebabTrd1Module.h"
68 //#include "AliEMCALRecPoint.h"
69 //#include "AliEMCALHistoUtilities.h"
70
71 ClassImp(AliEMCALGeometry)
72
73 // these initialisations are needed for a singleton
74 AliEMCALGeometry  *AliEMCALGeometry::fgGeom      = 0;
75 const Char_t*      AliEMCALGeometry::fgkDefaultGeometryName = "EMCAL_COMPLETEV1";
76 //
77 // Usage: 
78 //        You can create the AliEMCALGeometry object independently from anything.
79 //        You have to use just the correct name of geometry. If name is empty string the
80 //        default name of geometry will be used.
81 //         
82 //  AliEMCALGeometry* g = AliEMCALGeometry::GetInstance(name,title); // first time
83 //  ..
84 //  g = AliEMCALGeometry::GetInstance();                             // after first time
85 //
86 //  MC:   If you work with MC data you have to get geometry the next way: 
87 //  ==                                      =============================
88 //  AliRunLoader    *rl   = AliRunLoader::Instance();
89 //  AliEMCALGeometry *geom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
90 //  TGeoManager::Import("geometry.root");
91
92 AliEMCALGeometry::AliEMCALGeometry() 
93   : AliEMCALGeoUtils()
94
95   // default ctor only for internal usage (singleton)
96   // must be kept public for root persistency purposes, 
97   // but should never be called by the outside world    
98
99   AliDebug(2, "AliEMCALGeometry : default ctor ");
100 }
101 //______________________________________________________________________
102 AliEMCALGeometry::AliEMCALGeometry(const Text_t* name, const Text_t* title) 
103   : AliEMCALGeoUtils(name, title)
104 {
105   // ctor only for internal usage (singleton)
106   AliDebug(2, Form("AliEMCALGeometry(%s,%s) ", name,title));
107
108 }
109 //______________________________________________________________________
110 AliEMCALGeometry::AliEMCALGeometry(const AliEMCALGeometry& geom)
111   : AliEMCALGeoUtils(geom)
112 {
113   //copy ctor
114 }
115
116 //______________________________________________________________________
117 AliEMCALGeometry::~AliEMCALGeometry(void){
118     // dtor
119 }
120
121
122 //______________________________________________________________________
123 AliEMCALGeometry *  AliEMCALGeometry::GetInstance(){ 
124   // Returns the pointer of the unique instance
125   
126   AliEMCALGeometry * rv = static_cast<AliEMCALGeometry *>( fgGeom );
127   return rv; 
128 }
129
130 //______________________________________________________________________
131 AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
132                                                 const Text_t* title){
133     // Returns the pointer of the unique instance
134
135     AliEMCALGeometry * rv = 0; 
136     if ( fgGeom == 0 ) {
137       if ( strcmp(name,"") == 0 ) { // get default geometry
138          fgGeom = new AliEMCALGeometry(fgkDefaultGeometryName, title);
139       } else {
140          fgGeom = new AliEMCALGeometry(name, title);
141       }  // end if strcmp(name,"")
142       if ( AliEMCALEMCGeometry::fgInit ) rv = (AliEMCALGeometry * ) fgGeom;
143       else {
144          rv = 0; 
145          delete fgGeom; 
146          fgGeom = 0; 
147       } // end if fgInit
148     }else{
149         if ( strcmp(fgGeom->GetName(), name) != 0) {
150           printf("\ncurrent geometry is %s : ", fgGeom->GetName());
151           printf(" you cannot call %s ",name);  
152         }else{
153           rv = (AliEMCALGeometry *) fgGeom; 
154         } // end 
155     }  // end if fgGeom
156     return rv; 
157 }
158
159 //________________________________________________________________________________________________
160 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t distEff, Double_t &xr, Double_t &yr, Double_t &zr) const
161 {
162   // Jul 30, 2007 - taking into account position of shower max
163   // Look to see what the relative
164   // position inside a given cell is
165   // for a recpoint.
166   // In:
167   // absId   - cell is as in Geant,     0<= absId   < fNCells;
168   // e       - cluster energy
169   // OUT:
170   // xr,yr,zr - x,y,z coordinates of cell with absId inside SM 
171
172   // Shift index taking into account the difference between standard SM 
173   // and SM of half size in phi direction
174   const  Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
175   static Int_t nSupMod=0, nModule=-1, nIphi=-1, nIeta=-1, iphi=-1, ieta=-1;
176   static Int_t iphim=-1, ietam=-1;
177   static AliEMCALShishKebabTrd1Module *mod = 0;
178   static TVector2 v;
179   if(!CheckAbsCellId(absId)) return kFALSE;
180   
181   GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
182   GetModulePhiEtaIndexInSModule(nSupMod, nModule, iphim, ietam);
183   GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta); 
184   
185   //Get eta position. Careful with ALICE conventions (increase index decrease eta)      
186   if(nSupMod%2 == 0) {             
187     ietam = (fCentersOfCellsEtaDir.GetSize()/2-1)-ietam;// 47-ietam, revert the ordering on A side in order to keep convention.
188     if(nIeta == 0) nIeta = 1;
189     else           nIeta = 0;
190   }
191   mod = GetShishKebabModule(ietam);
192   mod ->GetPositionAtCenterCellLine(nIeta, distEff, v); 
193   xr = v.Y() - fParSM[0];
194   zr = v.X() - fParSM[2];
195   
196   //Get phi position. Careful with ALICE conventions (increase index increase phi)
197   Int_t iphi2 = iphi;
198   if(nSupMod<10) { 
199     if(nSupMod%2 != 0) 
200       iphi2 = (fCentersOfCellsPhiDir.GetSize()-1)-iphi;// 23-iphi, revert the ordering on C side in order to keep convention.
201     yr = fCentersOfCellsPhiDir.At(iphi2);
202     
203   } else {
204     if(nSupMod%2 != 0) 
205       iphi2 = (fCentersOfCellsPhiDir.GetSize()/2-1)-iphi;// 11-iphi, revert the ordering on C side in order to keep convention.
206     yr = fCentersOfCellsPhiDir.At(iphi2 + kphiIndexShift);
207   }
208   
209   AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
210   
211   return kTRUE;
212 }
213
214 //Not in use, comment for the moment
215 //________________________________________________________________________________________________
216 //Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Int_t maxAbsId, Double_t distEff, Double_t &xr, Double_t &yr, Double_t &zr) const
217 //{
218 //  // Jul 31, 2007 - taking into account position of shower max and apply coor2.
219 //  // Look to see what the relative
220 //  // position inside a given cell is
221 //  // for a recpoint.
222 //  // In:
223 //  // absId     - cell is as in Geant,     0<= absId   < fNCells;
224 //  // maxAbsId  - abs id of cell with highest energy
225 //  // e         - cluster energy
226 //  // OUT:
227 //  // xr,yr,zr - x,y,z coordinates of cell with absId inside SM 
228 //  
229 //  // Shift index taking into account the difference between standard SM 
230 //  // and SM of half size in phi direction
231 //  const  Int_t kphiIndexShift = fCentersOfCellsPhiDir.GetSize()/4; // Nov 22, 2006; was 6 for cas 2X2
232 //  static Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta;
233 //  static Int_t iphim, ietam;
234 //  static AliEMCALShishKebabTrd1Module *mod = 0;
235 //  static TVector2 v;
236 //
237 //  static Int_t nSupModM, nModuleM, nIphiM, nIetaM, iphiM, ietaM;
238 //  static Int_t iphimM, ietamM, maxAbsIdCopy=-1;
239 //  static AliEMCALShishKebabTrd1Module *modM = 0;
240 //  static Double_t distCorr;
241 //
242 //  if(!CheckAbsCellId(absId)) return kFALSE;
243 //  
244 //  GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
245 //  GetModulePhiEtaIndexInSModule(nSupMod, nModule, iphim, ietam);
246 //  GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi, ieta); 
247 //  
248 //  //Get eta position. Careful with ALICE conventions (increase index decrease eta)    
249 //  if(nSupMod%2 == 0) {             
250 //    ietam = (fCentersOfCellsEtaDir.GetSize()/2-1)-ietam;// 23-ietam, revert the ordering on A side in order to keep convention.
251 //    if(nIeta == 0) nIeta = 1;
252 //    else                 nIeta = 0;
253 //  }
254 //  
255 //  mod = GetShishKebabModule(ietam);
256 //  
257 //  if(absId != maxAbsId) {
258 //    distCorr = 0.;
259 //    if(maxAbsIdCopy != maxAbsId) {
260 //      GetCellIndex(maxAbsId, nSupModM, nModuleM, nIphiM, nIetaM);
261 //      GetModulePhiEtaIndexInSModule(nSupModM, nModuleM, iphimM, ietamM);
262 //      GetCellPhiEtaIndexInSModule(nSupModM,nModuleM,nIphiM,nIetaM, iphiM, ietaM); 
263 //      //Careful with ALICE conventions (increase index decrease eta)  
264 //      if(nSupModM%2 == 0) {             
265 //      ietamM = (fCentersOfCellsEtaDir.GetSize()/2-1)-ietamM;// 47-ietam, revert the ordering on A side in order to keep convention.
266 //      }
267 //      
268 //      modM = GetShishKebabModule(ietamM); // do I need this ?
269 //      maxAbsIdCopy = maxAbsId;
270 //    }
271 //    
272 //    if(ietamM !=0) {
273 //      distCorr = fEMCGeometry->GetEtaModuleSize()*(ietam-ietamM)/TMath::Tan(modM->GetTheta()); // Stay here
274 //      //printf(" distCorr %f | dist %f | ietam %i -> etamM %i\n", distCorr, dist, ietam, ietamM);  
275 //    }
276 //    // distEff += distCorr;
277 //  }
278 //  // Bad resolution in this case, strong bias vs phi
279 //  // distEff = 0.0; 
280 //  mod->GetPositionAtCenterCellLine(nIeta, distEff, v); // Stay here
281 //  xr = v.Y() - fParSM[0];
282 //  zr = v.X() - fParSM[2];
283 //  
284 //  //Get phi position. Careful with ALICE conventions (increase index increase phi)
285 //  Int_t iphi2 = iphi;
286 //  if(nSupMod<10) { 
287 //    if(nSupMod%2 != 0) 
288 //      iphi2 = (fCentersOfCellsPhiDir.GetSize()-1)-iphi;// 23-iphi, revert the ordering on C side in order to keep convention.
289 //    yr = fCentersOfCellsPhiDir.At(iphi2);
290 //    
291 //  } else {
292 //    if(nSupMod%2 != 0) 
293 //      iphi2 = (fCentersOfCellsPhiDir.GetSize()/2-1)-iphi;// 11-iphi, revert the ordering on C side in order to keep convention.
294 //    yr = fCentersOfCellsPhiDir.At(iphi2 + kphiIndexShift);
295 //  }
296 //  AliDebug(1,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
297 //  
298 //  return kTRUE;
299 //}
300 //
301
302 //
303 // == Shish-kebab cases ==
304 //
305
306 //
307 ////_________________________________________________________________________________
308 //void AliEMCALGeometry::GetGlobalEMCAL(const AliEMCALRecPoint *rp, TVector3 &vglob) const
309 //{
310 //  // Figure out the global numbering
311 //  // of a given supermodule from the
312 //  // local numbering for RecPoints
313 //
314 //  static TVector3 vloc;
315 //  static Int_t nSupMod, nModule, nIphi, nIeta;
316 //
317 //  const AliEMCALRecPoint *rpTmp = rp;
318 //  const AliEMCALRecPoint *rpEmc = rpTmp;
319 //
320 //  GetCellIndex(rpEmc->GetAbsId(0), nSupMod, nModule, nIphi, nIeta);
321 //  rpTmp->GetLocalPosition(vloc);
322 //  GetGlobal(vloc, vglob, nSupMod);
323 //}
324